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11.  AMTMCT  (AMmiinii  iOO  wordU 

An  investigation  was  performed  to  study  damage  development  in  fiber- reinforced  laminated 
composites  induced  by  transverse  concentrated  loads.  The  major  focus  of  the  study  was  to 
understand  fundamentally  the  damage  mechanics  in  terms  of  matrix  cracking  and  delami¬ 
nation  growth,  and  the  interaction  between  them  resulting  from  a  transverse  concentrated 
load.  The  study  was  focused  on  cross- ply  laminates  only,  and  the  load  was  introduced  quasi- 
statically  through  a  cylindrical  or  spherical  indenter.  Accordingly,  the  research  was  divided 
into  two  stages:  1)  damage  induced  by  a  cylindrical  indenter  and  2)  damage  induced  by  a 
spherical  indenter.  Analytical  models  consisting  of  a  contact  analysis  and  a  failure  analysis 
were  developed  for  analyzing  the  damage  initiation  and  growth  induced  by  both  loading 
conditions.  Experiments  were  also  performed  to  generate  data  needed  for  the  models  and 
to  verify  the  proposed  analyses.  The  predictions  based  on  the  models  agreed  very  well  with 
the  data. 
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Chapter  1 


Introduction 


It  is  well  known  that  fiber-reinforced  composites  are  vulnerable  to  transverse 
concentrated  loads  such  as  low- velocity  impact,  which  can  result  in  extensive  de¬ 
laminations  and  multiple  matrix  cracking  [1*34].  Such  internal  damage  can  cause 
significant  reduction  of  the  load  carrying  capacity  of  structmes  made  of  composites. 

Because  of  the  significance  of  the  problem,  numerous  experimented  and  analyt¬ 
ical  investigations  have  been  performed  to  study  the  damage  resulting  from  trans¬ 
verse  loads.  The  analytical  work  previously  developed  has  emphasized  estimating 
the  overall  delamination  sizes.  Wu  and  Springer  [4.  16]  proposed  a  non-dimensional 
empirical  expression  which  requires  several  parameters  to  be  determined  from  ap¬ 
propriate  experiments.  Gu  and  Sun  [21]  developed  a  model  for  estimating  the  im¬ 
pact  damage  size  in  sheet  molding  compound  (SMC)  composites.  The  applicability 
of  their  model  to  long  fiber  reinforced  laminates  is  questionable.  Choi  and  Chang 
[8,  32]  proposed  a  semi-empirical  model  based  on  failure  mechanisms  for  predicting 
the  delamination  size  and  locations  due  to  low- velocity  impact.  Although  the  model 
provides  reasonable  estimates  of  the  damage  size,  the  mechanics  of  damage  growth 
and  the  material  behavior  during  impact  cannot  be  obtained  from  the  model.  Finn 
and  Springer  [2]  developed  an  energy  criterion  for  predicting  delamination  sizes 
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and  locations  for  laminates  due  to  both  low-velocity  and  quasi-static  loading.  The 
model  needs  a  quantity  which  was  correlated  to  the  energy  per  unit  area  required 
to  delaminate  an  interface.  None  of  the  above  mentioned  models  considered  the 
detailed  fracturing  process  during  the  formation  of  the  damage. 

Several  investigators  have  indicated  based  on  their  experiments  that  for  some 
ply  orientations,  surface  matrix  cracking  could  be  produced  without  generating 
delaminations  by  carefully  adjusting  the  impact  velocity  or  the  amount  of  the  ap¬ 
plied  load.  However,  whenever  there  was  a  delamination,  there  were  matrix  cracks 
accompanied  with  the  delamination.  In  other  words,  matrix  cracking  and  delami¬ 
nation  were  detected  co-existing  in  most  test  conditions.  An  energy  threshold  also 
exists  below  which  no  damage  could  be  generated.  Experiments  have  demonstrated 
that  similar  results  obtained  from  low-velocity  impact  tests  could  be  produced  by 
quasi-static  transverse  loads  [1,  3,  7,  9]. 

Jones,  Joshi  and  Sun  [6,  7,  27]  have  suggested  that  matrix  cracking  is  the 
initial  impact  damage  mode.  Recently,  Choi  et  al.  [8,  32,  33,  34]  investigated 
impact  damage  using  a  cylindrical  indenter.  Because  the  indenter  produced  a  line 
load  during  impzict,  it  simplified  the  impact  damage  pattern  from  three  dimensions 
to  two  dimensions.  As  a  result,  the  damage  modes  could  be  detected  from  the  side  of 
the  specimens  by  the  nedced  eye  without  cutting  the  specimens.  Their  study  showed 
that  the  impact  energy  (velocity)  threshold  (  the  minimum  energy  required  to  cause 
damage)  is  associated  with  the  energy  required  to  initiate  matrix  cracks.  It  was 
concluded  that  matrix  cracking  is  the  initial  impact  damage  mode  and  delaminations 
are  generated  from  the  matrix  cracks. 

Baisically,  the  initial  matrix  cracks  have  been  classified  by  several  investigators 
[7,  8]  into  two  types:  the  shear  crack  and  bending  crack.  The  former  appears 
inside  the  laminate  and  is  located  a  distance  from  the  impacting  area,  but  the 
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latter  occurs  directly  beneath  the  transverse  load  in  the  outermost  layers  of  the 
laminate.  Previous  studies  [8,  32,  33,  34]  have  shown  that  the  extent  of  delamination 
strongly  depends  upon  the  location  of  the  initial  matrix  cracks.  Choi  et  al.  showed 
experimentally  from  their  line-loading  impact  study  that  the  shear  crack-induced 
damage  is  violent  and  catastrophic,  but  bending  crack-induced  damage  is  quite 
progressive.  It  is  apparent  that  there  is  a  strong  interaction  between  matrix  cracking 
and  delamination  growth  in  laminated  composites  resulting  from  transverse  loads 
[1-8].  In  contrast,  this  strong  interaction  does  not  appear  for  the  panel  due  to 
in-plane  loads. 

In  the  past,  matrix  cracking  and  delamination  have  been  frequently  treated  in 
the  literature  as  two  separable  damage  modes,  where  the  former  occurred  at  an  early 
stage  of  loading  and  was  related  to  the  initial  stage  of  failure,  while  the  latter  was 
associated  with  the  final  stage  of  failure  of  composites.  However,  in  transverse  load¬ 
ing,  matrix  cracking  once  developed  seems  to  affect  delamination  significantly.  This 
phenomenon  seems  to  set  apart  the  major  difference  in  deunage  growth  mechanisms 
between  transverse  concentrated  loadings  and  in-plane  loadings.  It  is  therefore  be¬ 
lieved  that  this  interaction  must  be  thoroughly  studied  in  order  to  fundamentally 
understemd  the  damage  mechanics.  Knowledge  of  the  damage  mechanics  is  criti¬ 
cally  important  to  material  scientists  for  improving  impact  resistance  of  composites 
and  to  design  engineers  for  selecting  and  designing  adequate  composite  structures 
which  may  be  subjected  to  transverse  loads. 

Limited  work  heis  been  performed  in  studying  the  damage  mechanics  and  the 
interaction  between  matrix  cracking  and  delamination  [38,  39,  46,  47,  48].  The 
models  of  Bostaph  and  Elber  [38]  and  Grady  et  al.  [39]  provide  an  estimate  of  the 
delamination  growth,  but  require  a  priori  knowledge  of  the  number  euid  locations 
of  delaminations.  Another  assumption  in  Bostaph  and  Elber’s  model  is  that  all 
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delaminations  grow  to  the  same  size.  Grady’s  model  applies  only  to  relatively  large 
delaminations.  However,  matrix  cracking  was  not  considered  in  their  reseeirch.  The 
fracture  modes  controlling  delamination  growth  of  laminated  composites  subjected 
to  transverse  concentrated  loads  have  not  yet  been  adequately  identified.  Further¬ 
more,  the  fracture  modes  contributing  to  the  delamination  propagation  could  be 
different  with  or  without  taking  into  account  the  effect  of  matrix  cracking.  Identi¬ 
fying  these  fracture  modes  will  be  important  for  designing  composite  material  and 
structural  systems  which  are  more  damage  resistroit  to  transverse  loads.  Sun  et  al. 
[46]  performed  a  2-D  linear  study  of  initiation  of  delamination  growth  induced  by 
a  surface  bending  crack  for  a  cross-ply  laminate.  They  found  that  Mode  I  domi¬ 
nates  the  initiation  of  delamination  growth.  The  very  recent  models  of  Martin  et 
al.  [47]  and  Salpekar  [48]  studied  the  effect  of  shear  crack  on  the  delzunination  for  a 
curved  panel  and  a  fiat  panel  respectively  by  fractiire  mechanics.  Their  2-D  mod¬ 
els  were  based  on  a  lineiir  analysis  and  only  initiation  of  the  delamination  growth 
was  considered.  Their  research  may  be  thought  to  be  parallel  to  the  2-D  model  of 
this  thesis  research  effort  [49-51].  Recently,  Liu  [55]  modeled  an  embedded  circular 
delamination  by  a  thin  soft  isotropic  layer  with  a  linear  finite  element  method  and 
showed  that  Mode  II  and  Mode  III  fractures  predominEuitly  control  the  growth  of 
the  embedded  delamination  due  to  a  point  load.  No  matrix  cracking  was  considered 
in  the  analysis. 

Accordingly,  the  objective  of  this  investigation  was  to  study  the  damage  me¬ 
chanics  of  laminated  composites  subjected  to  transverse  concentrated  loading  and 
to  fundamentally  understand  the  interaction  between  matrix  cracking  and  delami¬ 
nation  growth.  Both  flat  and  cylindrical  composites  were  considered.  To  simplify 
the  problem,  only  cross-ply  laminated  composite  panels  were  studied,  and  the  load 
was  applied  quasi-statically. 
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In  Chapter  2,  the  major  approach  for  achieving  the  objectives  of  the  investiga- 
;ion  is  briefly  outlined.  Chapter  3  describes  a  2-D  simplifled  study  both  analytically 
and  experimentally.  Simplified  damage  patterns  are  presented  and  damage  initia¬ 
tion  and  growth  is  discussed.  An  analytical  model,  consisting  of  a  stress  analysis, 
a  failure  analysis,  and  a  contact  analysis,  is  presented  in  detail  and  verified  both 
analytically  emd  experimentally.  The  interaction  between  matrix  cracking  and  de¬ 
lamination  is  discussed.  In  Chapter  4,  a  3-D  model  is  presented  to  study  the  effect 
of  initial  matrix  cracks  on  the  propagation  of  delamination  due  to  a  spherical  in- 
denter  loading.  The  model  zilso  consists  of  a  stress  analysis,  a  contact  analysis, 
and  a  failure  analysis.  The  contact  between  the  indenter  and  the  plate,  as  well  as 
between  the  delaminated  interfaces  is  modelled.  The  model  is  verified  extensively 
by  comparing  the  predictions  with  existing  data  and  experimentad  data.  Numerical 
examples  are  also  presented  for  a  toughened  composite  system  and  moderately  thick 
composites.  The  modes  contributing  to  the  delamination  propagation  are  identified 
and  discussed.  The  interfaces  for  the  2-D  and  3-D  computer  codes  are  presented  in 
Appendices  D  and  E. 


Chapter  2 


Method  of  Approach 


§2.1  The  Objectives 

This  research  was  performed  to  investigate  damage  development  in  fiber-reinforced 
lamainated  composites  induced  by  a  quasi-static  transverse  concentrated  load,  zis 
shown  in  Figure  2.1.  The  major  objectives  of  the  research  were  as  follows: 

1.  to  study  damage  mechanisms  and  mechauiics  of  fiber-reinforced  laminated  com¬ 
posites  due  to  quasi-impact  loading; 

2.  to  develop  adequate  models  for  predicting  quasi-impact  damage  in  the  materi¬ 
als’ 

3.  to  study  the  interaction  between  matrix  cracking  and  delsunination  growth. 

In  order  to  achieve  these  objectives,  this  research  was  ceirried  out  in  two  se¬ 
quential  studies  as  follows: 

1.  Damage  induced  by  a  cylindrical  indenter. 

2.  Damage  induced  by  a  spherical  indenter. 
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Figure  2.1  Description  of  the  problem.  A  laminated  composite  panel  subjected 
to  a  transverse  load  applied  by  a  spherically-nosed  indenter. 
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2.2.1  Damage  Induced  by  a  Cylindrical  Indenter 

A  2-D  analytical  model  was  developed  for  predicting  damage  in  laminated 
composites  resulting  from  a  quasi-static  line  load  applied  by  a  cylindrical  inden¬ 
ter.  Stresses  and  strains  in  composites  were  calculated  from  a  nonlinezir  finite 
element  code  based  on  generalized  plane  strsiin  conditions.  A  contact  model  was 
implemented  for  the  cracked/delaminated  systems  to  handle  contact  conditions  of 
cracked  interfaces.  A  modified  Hashin  matrix  failure  criterion  and  delamination  cri¬ 
terion  were  used  to  predict  the  initial  damage.  A  linear  fracture  mechanics  criterion 
was  adopted  for  predicting  the  damage  propagation.  Extensive  comparison  between 
the  predictions  and  experiments  was  made  during  the  research.  This  research  was 
expected  to  be  able  to  characterize  the  basic  failure  modes  and  mechanisms  of  lami¬ 
nated  composites  resulting  from  transverse  concentrated  loading  in  two  dimensions. 
The  focus  was  on  the  type  of  initizd  damage  mode  and  how  it  affected  the  damage 
growth.  The  stable  and  unstable  damage  growth  resulting  from  transverse  loading 
depending  on  the  location  of  the  initial  damage  were  emphasized. 

2.2.2  Damage  Induced  by  a  Spherical  Indenter 

A  3-D  analytical  model  was  also  developed  for  predicting  damage  in  laminated 
composites  resulting  from  a  quiisi-static  concentrated  load  applied  by  a  spherical 
indenter.  The  focus  was  on  the  effect  of  initial  matrix  cracking  on  the  delamination 
propagation.  The  model  consisted  of  a  stress  analysis  for  determining  the  stress 
distribution  inside  the  material  during  the  loading,  a  contact  einalysis  for  modeling 
the  contact  between  the  indenter  and  the  laminate  as  well  as  contact  between  the 
cracked  interfaces,  and  a  fEiilure  analysis  for  predicting  the  direction  and  extent  of 
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delamination  propagation.  In  the  stress  analysis,  finite  deformation  was  considered 
and  a  three-dimensionzJ  nonlinear  finite  element  was  used  to  provide  stress  and 
strain  information.  In  the  contact  analysis,  a  general  contact  model  was  proposed 
based  on  ein  augmented  Lagrangian  method  for  modeling  the  indenter-plate  contact 
and  cracked  interfaces  contact,  respectively.  Comparisons  were  also  made  between 
the  experiments  and  predictions  during  the  investigation. 


Chapter  3 


Damage  Induced  by  a  Cylindrical  Indenter 


§3.1  Statement  of  Problem 

In  this  problem,  a  flat  or  cylindrical  panel  made  of  fiber-reinforced  laminated 
composites  which  is  clamped  edong  two  parallel  edges  and  free  of  support  along  the 
others  was  considered.  A  transverse  concentrated  line  load  was  then  applied  at  the 
center  line  of  the  plate  parallel  to  the  two  clamped  edges,  as  shown  in  Figure  3.1. 
The  laminate  could  be  either  [0m/90„],  or  [90m/0n]»  layup,  where  m  emd  n  can  be 
arbitrary  integers.  It  wais  desired  to  obtain  the  following: 

1.  the  failure  and  the  response  of  the  structure  to  the  applied  load; 

2.  the  initiation  of  matrix  cracking  and  the  corresponding  magnitude  of  the  load; 

3.  delamination  initiation  and  propagation; 

4.  extent  of  delamination  m  a  function  of  the  applied  load. 

§3.2  Analysis 

The  analytical  model  proposed  for  describing  the  behavior  of  laminated  com¬ 
posite  panels  with  both  matrix  cracks  and  delaminations  consists  of  three  parts:  a 
stress  anedysis,  a  contact  analysis,  emd  a  failure  analysis.  The  stress  anedysis  was 
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developed  based  on  the  large  deformation  theory  for  analyzing  the  global  and  local 
deformations  of  the  sublaminates.  The  contact  analysis  used  the  Lagrange  multi¬ 
plier  method  for  modeling  the  matrix  crack  and  delamination  interface  condition 
during  loading.  The  failure  analysis  was  proposed  for  predicting  the  occurrence  of 
matrix  cracking  eind  for  modeling  delamination  initiation  and  propagation. 


3.2.1  Stress  Analysis 


As  the  local  deformations  of  the  laminate  could  be  substantial  due  to  the 
concentrated  load,  finite  deformation  theory  [70-74]  was  adopted  in  the  analysis. 
The  total  potential  energy  of  a  laminate  without  damage  under  the  given  load  caai 
be  described  as 


"pi*  r  f 

n  =  V  /  -  /  “T,  •  ui  °ds 


(3.1) 


where  rip/j,  is  the  total  number  of  layers,  is  the  cross-sectional  area  of  the 
mth  layer  in  the  reference  configuration,  and  °S  is  the  boundary  in  the  reference 
configuration  where  the  tractions  “T,  are  applied. 

The  components  of  the  second  Piola-Kirchhoff  stress  in  the  mth  layer  can  be 
expressed  by  the  constitutive  relation  as 


Qyn 

•^ij  —  ^ijkt  ^kl 


(3.2) 


where  .EJJ  axe  the  components  of  Green’s  strain  tensor.  are  the  orthotropic 

material  moduli  for  the  mth  layer. 

Based  on  the  virtual  work  principle,  the  following  equation  can  be  obtained 
from  Eq.  (3.1)  for  an  imcracked  laminate: 


2  f  S^SE^  “yn  -  /  “Ti  ■  8ui  ^ds  =  0 

m=i  J’S 


(3.3) 
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The  above  equation  is  valid  for  any  uncracked  laminated  composite  undergoing 
finite  deformations  in  both  two  and  three  dimensions. 

In  the  following,  the  formulation  is  given  for  a  laminated  composite  subjected 
to  a  line  load.  It  was  assiuned  that  the  width  of  the  laminate  was  considerably 
greater  than  its  thickness,  hence  the  problem  was  analyzed  two-dimensionally.  A 
two-dimensional  generalized  plane  strain  condition  in  the  xx—xz  plane  was  adopted. 
Therefore,  as  a  first  approximation,  the  displacement  field  of  the  l£Lminate  was 
assumed  in  the  form  [40-42]: 

Ui  =  Ui(Xi,l3) 

U2  =  ax2  (3.4) 

U3  =  U3(xi,X3) 

where  a  is  a  constant,  and  ui,U2,uz  are  the  displacements  in  the  x\  —  X2  —  xz 
coordinates,  respectively.  Consequently,  the  strain  E22  is  constant,  but  need  not  be 
zero.  However,  the  free  edge  effect  on  the  response  of  the  panel  is  neglected. 

Substituting  the  above  displacement  field  of  the  laminate  into  Eq.  (3.3)  provides 
a  nonlinear  equilibrium  equation  for  any  uncracked  laminate  subjected  to  a  line  load. 
However,  for  cracked  laminates,  Eq.  (3.3)  cannot  be  applied  directly  because  of  the 
presence  of  crack  interfaces  generated  inside  the  materials.  The  contact  condition 
of  the  interfacial  contact  during  loading  must  be  included  in  the  analysis. 

3.2.2  Contact  Analysis 

When  a  laminated  composite  is  subject  to  treinsverse  concentrated  loading,  the 
presence  of  matrix  crzwrks  and  delaminations  in  the  launinate  requires  that  the  con¬ 
dition  of  matrix  and  delamination  surfaces  during  deformation  be  specified  in  the 
analysis.  Because  of  material  mismatch,  relative  sliding  could  occur  along  the  in¬ 
terface  at  which  neighboring  plies  have  different  ply  orientation.  A  rigorous  contact 
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treatment  was  essential  for  accurate  modeling  of  the  interface  contact  conditions. 

In  order  to  prevent  the  contact  surfeices  from  overlapping,  an  impenetrability 
condition  must  be  specified  smd  satisfied  at  all  times  along  the  contact  interfaces. 
This  condition  requires  that  the  shortest  distamce  (defined  as  a  gap  g)  between 
two  contact  siirfaces  must  be  greater  them  or  equal  to  zero.  Mathematically,  the 
impenetrability  constraint  can  be  stated  as  ^  >  0.  Upon  contzict,  the  contact 
force  must  also  be  less  than  or  equal  to  zero  <  0).  Accordingly,  the  contact 
constraints  for  normal  contact  can  be  specified  as 

An  <0^1 

g{ui)  >  0  I  (3.5) 

An  =  0  J 

This  is  the  well-known  Kuhn-Tucker  conditions  for  normal  contact.  The  three  equa¬ 
tions  in  Eq.  (3.5)  reflect,  respectively,  the  compressive  normal  traction  constraint, 
the  impenetrability  constraint,  and  the  requirement  that  the  pressure  is  nonzero 
only  when  ^  =  0. 

Accordingly,  Eq.  (3.5)  must  be  considered  together  with  Eq.  (3.3)  in  order  to 
analyze  composites  containing  interned  cracks  and  delaminations.  Several  methods 
exist  for  implementing  the  contact  constraints  [62-63].  One  of  the  most  popular 
choices  is  the  Lagrange  multiplier  method  [63].  The  Lagrange  multiplier  method 
enforces  the  exact  constraints,  and  it  was  adopted  in  the  present  two-dimensioned 
study. 

By  imposing  the  constraint  conditions  into  the  total  potential  energy  using  the 
Lagrange  multiplier  technique,  the  total  potential  energy  given  in  Eq.  (3.1)  can  be 
modified  as 


n( 


Ui,AN)=  V  /  f  ^Ti-urdS+  f  XsgMdr  (3.6) 

J»s  Jr 
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where  is  the  Lagrange  multiplier  and  F  is  the  sum  of  matrix  crack  and  delami¬ 
nated  surfaces  in  contact. 

The  variational  equations,  corresponding  to  the  Lagrange  multiplier  formula¬ 
tion,  result  in 


V  f  ST,  SET,  “dil  -  f  °Ti-  Sui  ^ds  +  /  XsSgdT  =  0,  ] 

Jt  I 

j  SXfjg  dT  =  0 

Jr 

Or  in  a  more  concise  form,  the  above  equations  may  be  rewritten  as 

=  0, 

J  dT  =  0 


(3.7) 


(3.8) 


where  ^11  is  the  virtual  work  for  an  unconstrained  system,  zind  SXfj  and  Sg  are  the 
variation  of  A;v  and  g,  respectively. 

In  order  to  solve  the  above  equations,  the  gap  fu.iction  g  needs  to  be  defined 
specificeilly  for  any  contact  interface.  For  the  purpose  of  generalization,  ziny  one  of 
the  contact  bodies  czin  be  defined  as  “contactor”  emd  the  other  can  be  referred  to  as 
“target”.  Any  points  belonging  to  the  contact  region  of  the  contactor  will  be  called 
“slave  points”  and  those  on  the  surface  of  the  target  will  be  referred  to  as  “meister 
points”.  In  the  following,  m  and  s  indicate  points  on  the  master  surface  and  the 
slave  surface,  respectively. 

Therefore,  for  a  given  pair  of  potenti2il  contact  points  2dong  the  contact  surfaces, 
aj,  on  the  contactor’s  surface  and  Xc  on  the  target’s  surface,  the  gap  function  is 
defined  as 


ls(«i)l  =  ll®»  -  ®ci!  =  min  ||X,  -  Xml 


(3.9) 


where  Xm  cam  be  any  point  on  the  target  surface.  Expressed  in  words,  Xc  (contact 
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Figure  3.2  Local  coordinate  systems  of  contact  surfaces  in  the  two-dimensional 


study. 
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point)  is  the  closest  point  projection  of  x,  onto  the  target  snrf2w:e,  as  shown  in 
Figure  3.3. 


It  is  necessary  to  give  an  accurate  description  of  the  local  geometry  on  the 
contact  surfaces  during  the  deformation.  Here,  a  local  coordinate  system  needs  to 
be  defined  such  that  it  is  suitable  for  a  definition  of  the  contact  constraints.  This 
local  coordinate  system  is  fixed  along  the  target  surface,  as  shown  in  Figure  3.3. 
Before  the  local  coordinate  system  can  be  established,  it  is  necessary  to  parametrize 
the  contact  surface  locally.  In  terms  of  local  surface  parameters  the  following  basis 
for  any  master  point  on  the  target  surface  is  defined: 


C.{«)  =  (or 

C2(0  =  (0,1,0)^  or  (0,-l,0)^ 

g  /tx  _  fi(e\  —  ^  ^2(0 

l|ei(f)  X  e,(01l  ) 


(3.10) 


where  61,62,  and  63  are  vectors  in  3-D  Cartesian  space.  The  sign  of  62  depends 
on  the  definition  of  the  target  surface  such  that  the  normal  Tl  is  outward  on  the 
target  surface  (Figure  3.3). 

Accordingly,  the  gap  function  g  can  be  written  as 


g  =  n-d  =  n-{x,-Xc)  (3.11) 

This  definition  is  true  for  both  3-D  and  2-D  contact  gap  definition.  In  terms  of  local 
coordinate  system  the  gap  g  can  be  written  as  follows: 

g  =  {x>-Xc)-n{^c)  (3.12) 

In  order  to  solve  Eq.  (3.8),  a  nonlinear  finite  element  method  was  developed. 
The  numerical  formulation  will  be  presented  in  section  3.3. 
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g  <  0  :  overlapping 
g  =  0  :  contact 

g  >  0  :  no  contact 


Figure  3.3  The  definition  of  gap  function  in  two-dimensional  space. 
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The  primary  feulure  modes  of  interest  were  matrix  cracking  and  delamination, 
hence  fiber  breakage  was  not  considered.  Failure  criteria  were  proposed  for  pre¬ 
dicting  the  occurrence  of  the  initial  failure.  A  fracture  analysis  b^lsed  on  the  linear 
elastic  fracture  mechanics  theory  was  also  adopted  for  modeling  the  crack  propaga¬ 
tion  once  the  damage  initiated. 


Prediction  of  Initial  Failure 

In  order  to  predict  initial  damage,  a  matrix  failure  criterion  and  a  delzimination 
failure  criterion  [32,  44]  were  adopted.  The  matrix  failure  criterion  can  be  described 
as 


(3.13) 


The  occurrence  of  matrix  cracking  is  predicted  when  the  value  of  is  equal  to  or 
greater  than  unity. 

The  delamination  failure  criterion  can  be  expressed  as 

(3.14) 


» *  (^)  ■ 


The  occurrence  of  delamination  at  the  interface  was  predicted  when  the  value  of 
was  equed  to  or  greater  than  unity.  In  Eq.  (3.13)  and  Eq.  (3.14),  Si  is  the  in- situ 
interlaminar  shear  strength  within  the  laminate  under  consideration,  and  Yt  is  the 
in-situ  ply  transverse  tensile  strength  [32].  <7y,  and  <7„  are  the  interlaaninax  shear 
stresses,  and  and  Cyy  are  the  out-of-plane  normal  stress  and  in-plane  transverse 
normal  stress  respectively  within  the  laminate  under  consideration. 

Whenever  the  combined  state  of  stresses  satisfied  either  one  of  the  criteria, 
initial  failure  was  predicted.  The  corresponding  failure  criterion  indicated  the  initizd 
mode  of  failure.  Once  the  initial  failure  was  predicted,  fracture  emzdysis  was  applied 
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to  simulate  the  growth  of  the  local  damage  as  the  applied  load  continued  to  increase. 
Modeling  of  Crack  Propagation 

In  order  to  simulate  crack  propagation,  a  small  initial  crack  was  introduced 
immediately  after  the  occurrence  of  the  initial  failure,  depending  upon  the  type  of 
the  failure  mode.  For  matrix  cracking,  a  vertical  crack  was  immediately  generated 
in  the  failed  ply  at  the  location  where  matrix  cracking  was  predicted.  The  size  of 
the  crack  was  assumed  to  be  equal  to  the  thickness  of  the  cracked  ply  group  and 
fully  extended  to  the  interfaces  of  the  neighboring  plies,  which  have  different  ply 
orientations.  However,  if  delamination  was  predicted  as  the  initial  failure  mode, 
an  initial  delamination  length  about  one  ply  thickness  would  be  introduced  at  the 
predicted  interface.  The  stresses  and  deformations  of  the  laminate  would  then  be 
recalculated  at  the  same  load,  and  fracture  analysis  would  be  applied  to  determine 
the  growth  of  the  initial  failure  and  to  model  the  propagation  of  the  damage. 

The  mixed  mode  fracture  criterion  that  was  adopted  to  predict  the  initiation 
of  crack  propagation  can  be  expressed  as  [42,  54,  53]: 


(3.15) 


where  Gic  and  Gjic  are  the  criticed  strain  energy  release  rates  corresponding  to 
Mode  I  and  Mode  II  fracture,  respectively.  It  was  assumed  that  the  values  of  Gjc 
and  Giic  did  not  change  with  delamination  length,  a  =  1  and  0  =  1  were  selected 
for  this  study  because  they  were  found  to  provide  the  best  fit  to  the  previous 
experiments  [42]. 

The  onset  of  crack  growth  was  predicted  when  the  value  of  was  equal  to  or 
greater  than  unity  {Ed  >  1).  The  strain  energy  release  rates  G/  and  G//  for  Mode 
I  and  Mode  II  fractures,  respectively,  could  be  expressed  based  on  the  linear  eleistic 
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fracture  mechanics  as 

(  1 

Gj  =  lim  \ 

Aa— .0 

[^L 

[u;|‘(a  -1-  Aa)  -  u„  (a  -f-  Aa)]  <t„„  ds  > 

(3.16.a) 

f  j  ^Aa 

r  1 

Gii  =  lim  ^ 
Aa— .o 

[^lo 

[u;  (a  +  Aa)  -  u,  (a  -}-  Aa)]  (r,„  ds  > 

(3.16.6) 

where  Aa  is  the  crack  extension,  n 

and  s  denote  normal  zmd  sliding  direction,  (7„n 

and  (7n«  are  the  stresses  at  the  crack  tip  associated  with  a  delamination  length  of 
a,  and  ujJ",  it^,  and  u~,  uj  are  the  displacements  of  the  upper  and  lower  surfaces 
of  the  delamination  associated  with  a  delamination  length  of  a  +  Ac  respectively. 

§3.3  FINITE  ELEMENT  MODELING 

A  nonlinear  finite  element  code  was  developed  for  implementing  the  proposed 
model.  The  laminate  was  discretized  by  4-node  bilinear  elements  into  Uti  element 
domains  with  n„p  nodal  points  using  isoparametric  4-node  bilinear  rectangular  el¬ 
ements  as  shown. 

The  displacement  field  is  assumed  to  be  in  the  form 

U  =  5^iV.Ua  (3.17) 

asl 

where  ne„  is  the  number  of  nodes  for  one  element,  Na  is  the  displacement  shape 
(interpolation)  function  associated  with  the  local  node  number  a,  and  Ua  is  the 
nodal  displacement  vector  at  the  local  node  number  a.  The  displacement  shape 
functions  in  the  parent  domain  with  coordinates  p  and  q  are  expressed  as 

^a(p,9)=  J(l+PP«)(l  +  99a)  (3.18) 

where  a  =  1  —  4  emd  there  is  no  sununation  on  a. 

Now  that  both  contact  bodies  are  discretized  by  finite  elements,  any  one  of 
the  contact  bodies  can  be  assigned  as  the  contactor  and  the  other  will  be  treated 
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as  the  target.  The  nodal  points  on  the  contact  region  of  the  contactor  are  ceJled 
“slave  nodes,”  and  those  on  the  surface  of  the  tzirget  are  called  “master  nodes”.  The 
elements  containing  the  master  nodes  are  called  “master  elements,”  and  the  those 
containing  the  slave  nodes  are  called  “slave  elements”.  The  master-slave  definition 
is  shown  in  Figure.  3.4. 

In  the  current  formulation,  the  Lagrange  multiplier  distribution  is  assumed  to 
be  constant  for  a  slave  node  associated  with  a  master  element  along  a  crack  or 
delamination  interface.  Thus,  the  extra  term  in  the  total  potentizil  energy  due  to 
the  presence  of  Lagrange  multipliers  can  be  expressed  as 

/  Xrsf  gdT  =  f  {\N)igidri  (3.19) 

1^1  •'r. 

where  n,  is  the  total  number  of  slave  nodes  and  therefore,  the  number  of  contact 
elements  on  the  contact  surfaces.  Because  of  the  constant  contact  pressure  assump¬ 
tion  for  each  contact  element  and  the  unique  g  for  each  contact  element,  the  above 
equation  can  be  written  as 


f  XNgdT^TiXNygih  (3.20) 

where  li  is  the  length  of  the  i-th  contact  element.  In  the  discrete  space,  the  sum 
of  the  contact  pressure  {Xj^l)-  must  be  carried  by  slave  node  s,  which  could  be 
imagined  8is  a  nodal  force  with  a  magnitude  of  (A;v0i-  Thus,  the  extra  term  in  the 
total  potential  energy  due  to  the  presence  of  Lagrange  multipliers  can  be  expressed 
as 


/  X^gdr  =  y^  Xigi 


(3.21) 


where  the  (=  (A^vOi)  ^9-  (3-21)  is  interpreted  as  the  Lagrange  multi¬ 

plier  for  a  slave  node  along  some  matrix  crack  or  delamination  interface. 
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Master  Surface 


Figure  3.4  Description  of  master  £ind  slave  surfaces  for  the  contact  problem 


two-dimensioned  space. 
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3.3.1  Vziriational  Formulation  of  the  Discrete  Problem 

Accordingly,  Eq.  (3.6)  can  be  written  in  its  discrete  form  as 


n(u,A)  =  n(u)  +  A'^g 


(3.22) 


where  u  is  the  vector  of  nodal  point  displacements  in  discrete  space,  g  the  vector  of 
nodal  gaps,  and  A  the  vector  of  nodal  contact  forces.  For  simplicity,  the  subscript  N 
is  omitted  for  A.  In  addition,  n(u)  designates  the  total  potential  energy  associated 
with  the  bodies  in  contact.  The  discrete  form  of  the  virtual  work  principle  results 


in 


dU  dU  9TI  dU  t  c  x 

^  =  <  a;  +  ^  + «  "  =  “ 


(3.23) 


Therefore,  the  nonlinear  equilibrium  equations  for  the  discrete  system  are  ob¬ 
tained  as 


g  =  o  J 


(3.24) 


In  the  framework  of  the  finite  element  method,  the  term  can  be  replaced  by 


an  _ 

^  -  Pint  -  R.xt 


(3.25) 


where  Pint  is  the  vector  of  internal  nodal  forces  resulting  from  the  element  stresses 
(corresponding  to  the  strain  energy  part  in  the  total  poential  energy  H),  and  R^xt 
is  the  vector  of  nodal  forces  due  to  external  loading.  Term  A^Jj  can  be  replaced 
by 


‘'I-'- 


(3.26) 


where  the  F*  is  the  vector  of  contact  nodal  forces  due  to  the  contribution  of  contact 
nodes. 
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Combining  Eq.  (3.24),  Eq.  (3.25),  and  Eq.  (3.26)  results  in 

Pint  "I"  Pc  —  ^ 

g  =  0 


(3.27) 


The  first  equation  represents  a  weak  form  of  equilibrium  equation.  The  second 
indicates  that  the  constraints  are  exactly  enforced  here.  It  is  shown  that  A  represents 
the  vector  of  contact  forces,  which  are  distributed  according  to  the  matrix  in 
the  associated  nodes. 

In  general,  due  to  the  finite  deformation  and  the  contact  constraints,  the  sys¬ 
tem  of  equations  is  nonlinear.  NumericzJ  techniques  such  as  the  Newton-Raphson 
method  have  been  widely  used  to  solve  these  nonlineeir  equations  [70-74]. 

By  applying  the  Newton-Raphson  technique,  the  following  linearized  equations 
are  obtained: 


^[Pinl  +  Fcl<fu  +  ^[Pi«  +  F.HA  =  R„,  -  (Pin.  +  F.), 


5u 


(3.28) 


du  =  -g 


Recalling  that  Pjnt  is  the  fimction  of  u  only  and  Fc  is  the  function  of  u  and 
A,  the  above  equation  can  be  written  as 


l^(Pin.)  +  |;(Fc)|rfu  +  l|f(F.)lrfA  =  Rnx.  -  (Pin.  +  F.), 


5u 


du  =  -g 


The  above  equations  can  be  written  in  matrix  form: 


(3.29) 


(K  Kc)du  -H  HdA  =  R.xt  -  (Pint  +  Fc), 

H^du  =  — g 


1 


(3.30) 


where  is  the  transpose  of  H. 
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In  this  matrix  formulation,  K  is  the  classical  tangent  stiffness  matrix  which 
is  obtained  without  considering  the  contact  between  crack  interfaces,  and  is  the 
standard  practice  in  finite  element  formulation  [70-74]  amd  therefore  omitted  here. 
Fc  is  the  contact  residual  force  vector  and  defined  as 


Fc 


(3.31) 


Ko  is  the  contribution  of  the  contact  conditions  to  the  tangent  stiffness  matrix, 
which  guarantees  the  best  approximation  to  the  actuad  response.  It  is  defined  as 


Kc  = 


aFe 

du 


(3.32) 


H  is  a  vector  that  sets  up  the  distribution  of  contact  forces  due  to  the  contact 
and  is  defined  as 


(3.33) 


In  the  following  section,  it  is  convenient  to  refer  the  equations  to  an  element 
coming  in  contact  with  one  node  on  the  other  body,  forming  one  contact  element. 
Such  an  element  provides  one  gap  g  and  one  Lagrange  multiplier  A  eind  has  N  nodal 
points.  Subsequently,  capital  letter  indices  are  used  here  to  identify  quantities  which 
are  related  to  the  element  nodal  points.  Moreover,  for  repeated  indices  the  summa¬ 
tion  convention  of  tensor  calculus  holds.  The  corresponding  form  of  Eq.  (3.30)  for 
one  contact  element  is 


I  \  ^  r  11  I  - 


(PliJt+F-), 


dg 

du^ 


dn^  =  -g 


(3.34) 


where  the  contact  residual  force  of  the  element  is 
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The  two  terms  in  parentheses  form  the  tangent  stiffness  matrix  of  the  element. 
The  first  term  is  the  classical  tangent  stiffness  matrix  for  imdamaged  composites. 
The  derivation  of  the  matrix  is  available  in  the  literature  related  to  finite  element 
methods  [70-74].  The  second  term  adding  to  the  contribution  due  to  contact  has 
the  following  form: 


k  =  ik  = 


(3.36) 


auN  auN'auM 

The  correct  derivation  of  contact  element  stiffness  matrix  kc  is  important,  as 
it  renders  the  quadratic  convergence  of  Newton’s  method.  The  coefficient  matrix 
before  the  dX  forms  the  contact  force  distribution  of  the  associated  nodes  of  the 
contact  element.  The  derivations  of  the  contact  stiffness  matrix  kc  and  the  residual 
force  t  are  described  in  Appendix  B. 


3.3.2  Finite  Element  Modeling  of  Failure 

Once  the  strain  and  stress  distributions  in  the  laminate  at  each  load  increment 
were  calculated,  failure  criteria  were  applied  to  determine  the  locations  of  the  initial 
damage  and  the  applied  load  corresponding  to  the  initiad  failure.  When  intraply  or 
interply  cracking  was  predicted,  it  was  necessary  to  insert  a  crack  in  remodeling  the 
finite  element  grids  within  and  near  the  cracked  area.  However,  because  the  location 
of  the  crack  weis  unknown  prior  to  failure  and  also  strongly  dependent  upon  the  ply 
orientation,  it  became  very  difficult  to  efficiently  insert  ein  intraply  or  interface  crack 
in  the  existing  mesh  and  effectively  remesh  the  finite  element  grids  while  the  ply 
orientation  and  the  geometry  of  the  laminate  were  considered  as  variables. 

Therefore,  an  efficient  and  effective  technique  weis  developed  and  implemented 
in  the  finite  element  mesh  generator.  The  technique  used  the  original  finite  element 
grid  system,  amd  preserved  the  sparseness  of  the  tangent  stiffness  matrix  after  failure 
occurred.  Hence  the  computational  time  did  not  increase  substantially.  When  a 
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crack  was  predicted,  it  was  assumed  that  the  location  of  the  crack  could  be  assigned 
8LS  an  approximation  to  the  nearest  element  boundary.  Accordingly,  a  matrix  crack 
with  a  length  of  a  ply  group  thickness  would  be  generated  along  the  existing  element 
interfaces  of  the  failed  ply  group.  Similarly,  for  a  delamination  crack,  the  interfacial 
crack  would  be  introduced  at  the  boundary  of  the  elements  sdong  the  predetermined 
ply  interface. 

Based  on  this  approach,  each  node  was  systematically  assigned  four  node  num¬ 
bers  (rather  than  one)  in  the  global  nodal  numbering  system.  If  no  damage  occurred, 
the  four  nodes  were  assigned  the  same  equation  numbers.  Once  a  crack  was  gen¬ 
erated,  depending  upon  the  type  of  crack  and  the  location  of  the  node,  each  node 
along  the  cracked  surfaces  could  be  generated  with  its  own  nodal  number,  up  to 
four  individual  nodes.  Accordingly,  additional  crack  surfaces  could  be  generated  by 
splitting  the  nodes  at  the  desired  locations.  Hence,  no  changes  in  the  mesh  were 
required  for  the  post-failure  analysis.  Details  of  the  mesh  remodeling  procedure  are 
given  in  Appendix  A. 

Once  a  crack  was  introduced  into  the  composite,  the  surface  condition  of  the 
crack  during  the  subsequent  loading  had  to  be  modeled  properly  to  prevent  over¬ 
lapping  of  the  crack  surfaces.  It  was  assumed  that  once  matrix  cracking  occurred, 
the  crack  would  fully  extend  to  the  neighboring  interfaces  of  plies  with  ply  orien¬ 
tations  differing  from  the  cracked  ply  group  and  would  branch  into  a  delamination 
crack  along  each  interface.  In  order  to  apply  fracture  mechanics  for  determining  the 
growth  of  the  matrix  cracking-induced  delamination,  ztn  initial  deleunination  crack 
along  the  interface  perpendicular  to  the  matrix  crack  was  generated  at  the  tips  of 
the  existing  matrix  crack.  Thus,  the  strmn  energy  releeise  rates  at  the  delamination 
tips  along  the  interface  could  be  determined  based  on  the  crack  closure  technique 
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(Eqs.  3.16)  [48]  as  follows: 

“  ii-o  { 

*  ii“o  { lij 

where  M  and  N  are  the  nodes  at  the  crack  tip  before  the  crack  advances  an  amo\int 
of  Aa,  FnM  and  F,m  are  tlie  nodal  forces  that  keep  the  nodes  M  and  N  together, 
and  Aa  is  chosen  to  be  the  length  of  the  element  sJiead  of  the  crack  tip  along  the 
delamination  interface. 

It  was  assumed  that  the  matrix  cracking-induced  delamination  crack  would 
propagate  if  the  calculated  strain  energy  release  rates  Gj  and  Gu  satisfied  the 
crack  growth  criterion  of  Eq.  (3.15).  Subsequently,  the  delamination  crack  would 
axlvsmce  by  one  element  size,  and  the  mesh  with  additional  crack  surfaces  would 
require  remodeling  following  the  procedures  given  in  Appendix  A.  Accordingly, 
Eq.  (3.27)  needed  to  be  resolved  with  additional  crack  surfaces.  The  above  numerical 
procedures  would  need  to  be  continued  until  either  the  composite  could  not  sustain 
any  additional  load  or  the  calculated  results  would  be  beyond  the  scope  of  the 
present  interest. 


§3.4  Verification 

In  order  to  verify  the  model,  numericiil  calculations  based  on  the  model  were 
compared  with  the  available  analytictil  and  experimental  results.  Experiments  were 
also  performed  during  the  investigation  to  further  validate  the  model. 

3.4.1  Comparison  with  Existing  Results 
A  DCB  Specimen  under  a  Line  Load 

Numerical  calculations  were  performed  to  determine  the  strain  energy  release 
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rates  of  a  double  cantilever  beam  subjected  to  a  constant  tip  displacement  under 
three  different  loading  conditions.  The  results  of  the  calculations  are  presented  in 
Figure  3.5. 

Figure  3.5  shows  the  strain  energy  release  rate  calculations  for  Mode  I,  Mode 
II,  and  mixed  mode  fracture,  as  a  function  of  the  crack  length  under  a  constant 
applied  displacement  at  the  crack  tip.  The  solid  lines  represent  the  results  of  the 
beam  theory  and  the  symbols  are  the  predictions  based  on  the  present  analysis.  The 
predictions  agreed  very  well  with  the  classical  beam  theory  [45].  For  mixed  mode 
fracture,  Gj  and  Gu  were  calculated  for  the  model  and  are  presented  in  Figure  3.5, 
in  addition  to  the  total  strain  energy  release  rate  G  =  Gi  -¥  Gu. 

A  Unidirectional  Composite  under  Three-Point  Bending 

Calculations  were  made  to  compare  the  prediction  with  the  test  data  (Fig¬ 
ure  3.6)  measured  from  experiments  in  which  a  unidirectional  composite  beam  con¬ 
taining  a  central  delamination  was  tested  imder  three-point  bending  (by  Maikuma 
et  al.  [43]).  The  Strain  energy  release  rates  with  respect  to  the  crack  length  were 
measured.  The  open  symbols  represent  the  test  data,  and  the  solid  line  represents 
the  predictions  from  the  model.  The  predictions  agreed  with  the  test  data  reason¬ 
ably  well. 

[90,„/0„/90m]  Specimens  with  a  Surface  Crack 

Comparisons  were  also  made  between  the  predictions  beised  on  the  model  and 
the  test  data  measured  by  Sun  et  al.  [46]  on  a  simply  supported  [9O5/O5/9O5]  com¬ 
posite  beam  containing  an  interface  delamination  induced  by  a  surface  matrix  crack 
(bending  crack)  located  at  the  center  of  the  bottom  ply.  The  test  was  performed 
to  evaluate  the  e  ffect  of  a  bending  crack  on  the  strain  energy  release  rates  at  the 
tip  of  a  delamination.  Figure  3.7  presents  the  comparison  of  strain  energy  release 
rates  between  the  calculated  rates  based  on  the  model  and  those  given  in  [46].  A 


CRACK  TIP  LOCATION  (i/L) 


Figure  3.5  Strain  energy  release  rates  of  an  isotropic  beam  as  a  function  of  the 
crack  length.  Comparison  is  between  the  predictions  based  on  the 
present  model  and  the  calculations  based  on  the  classical  beam  theory 
[45].  a)  Mode  I;  b)  Mode  II;  and  c)  Mixed  Mode. 


ENERGY  RELEASE  RATE  G/P^  (l.OE-4/Nm) 
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Figure  3.6  Comparison  of  the  strain  energy  relesise  rate  of  an  AS4/2220-3  [O24] 
composite  beam  between  the  data  taken  from  [43]  and  the  predictions 
based  on  the  present  model. 
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fairly  good  agreement  was  found.  It  was  noticed  that  Mode  I  was  the  dominant 
mode  for  the  crack  growth.  Figure  3.8  shows  the  relationship  between  the  critical 
load  associated  with  the  onset  of  the  delamination  growth  and  the  delamination 
length  for  the  launinate  under  a  three  point  bending  test.  The  rectangular  symbol 
represents  the  test  data.  The  solid  line  is  the  prediction.  The  prediction  was  based 
on  the  assumption  that  the  critical  strain  energy  release  rates  Gjc  and  Giic  were 
constant  for  the  given  material  and  independent  of  the  delamination  length.  Again, 
the  predictions  agreed  with  the  data  very  well. 

The  nximericaJ  simiilations  of  the  deformed  configurations  of  a  [9O8/O4],  com¬ 
posite  specimen  with  a  bending  crack  as  a  function  of  the  applied  load  are  pre¬ 
sented  in  Figure  3.9.  The  calculated  load  amd  the  interface  delamination  versus 
the  displacement  of  the  load  head  were  also  calculated  and  are  presented  in  Fig¬ 
ure  3.10.  Clearly,  the  interfacial  delamination  induced  by  the  bending  crack  grew 
gradually  with  the  increasing  applied  load  or  displacement,  indicating  a  very  stable 
crack  growth.  The  strain  energy  releeise  rates  were  calculated  and  axe  presented 
in  Figure  3.11,  which  shows  again  that  Mode  I  wzis  the  dominemt  mode  during 
delamination  growth. 


ENERGY  RELEASE  RATES  (G/P>  )  (1/Nm) 
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0.0  0.2  0.4  0.6 

CRACK  TIP  LOCATION  (a/L) 


Figure  3.7  Strain  energy  release  rate  of  an  AS4/3401-6  [9O5/O5/9O5]  composite 
beam  as  a  function  of  the  delamination  length  induced  by  a  surface 
matrix  crack.  Comparison  is  between  the  predictions  based  on  the 
present  model  and  the  calculations  given  in  [46]. 


CRITICAL  LOAD  Per  (lbs) 
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CRACK  TIP  LOCATION  (a/L) 


Figure  3.8  Comparison  between  the  measured  [46]  and  the  predicted  critical  loads 
of  an  AS4/3401-6  [90s/0s/90s]  composite  beam  as  a  function  of  the 
delamination  length  induced  by  a  surface  matrix  crack. 
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P  =  0.0  lbs 


1300,976  [908,04]s 


Figure  3.9  Numerical  simulations  of  the  deformed  configurations  of  a  T300/976 


[9O8/O4],  composite  beam  containing  an  existing  matrix  crack  and 
subjected  to  a  transverse  concentrated  line  load  at  various  loading 
stages. 
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Figure  3.10  The  calculated  response  of  a  T300/976  [9O8/O4],  composite  beam  con¬ 
taining  a  bending  craek  and  subjected  to  a  transverse  concentrated  line 
load.  (Above);  The  calculated  applied  load  as  a  function  of  the  load 


head  displacement.  (Below):  The  extent  of  the  delamination  induced 


by  the  matrix  crack  as  a  function  of  the  load  head  displacement. 
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DISPLACEMENT  A  (in) 


Figure  3.11  The  ratios  of  the  calciilated  strain  energy  release  rates  corresponding 
to  the  delamination  growth  as  a  function  of  the  load  head  displace¬ 


ment. 
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Table  3-1 

Material  properties  used  in  the  calculations  for  comparisons. 


Aluminum 

T300/976 

AS4/3501-6 

AS4/2220-3 

E^,  (msi) 

10.0 

17.6 

17.4 

17.84 

Ey,  (msi) 

10.0 

1.41 

1.43 

1.61 

Ez,  (msi) 

10.0 

1.41 

1.43 

1.61 

t'xy 

0.3 

0.29 

0.30 

0.29 

Vxz 

0.3 

0.29 

0.30 

0.29 

Vyz 

0.3 

0.40 

0.30 

0.40 

Gxy,  (msi) 

3.846 

0.81 

0.76 

0.908 

Gzz,  (msi) 

3.846 

0.81 

0.76 

0.908 

Gyz,  (msi) 

3.846 

0.50 

0.76 

0.623 

Gi^,  (Ibf/in) 

n/a 

0.50 

0.47 

n/a 

Guc,  (Ibf/in) 

n/a 

1.80 

1.88 

2.80 

For  the  materials  T300/976,  AS4/3501-6  and  AS4/2220-3,  it  is  necessary  to 
specify  the  direction  in  which  the  fibers  of  each  lamina  axe  aligned.  The  subscript 
X  denotes  the  fiber  or  Odeg  direction.  The  subscript  y  denotes  the  direction  per¬ 
pendicular  to  the  Odeg  direction  (transverse  or  90 deg  direction).  The  subscript  z 
denotes  the  out-of-plane  direction. 
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3.4.2  Comparison  with  the  Current  Test  Data 

Tests  were  also  performed  during  the  investigation,  and  the  results  of  the  test 
data  were  utilized  to  substantiate  the  model.  The  major  focus  of  the  tests  W2is 
concentrated  on  damage  in  the  cross-ply  Izuninates  [0m/90„],  contedning  inner  90° 
plies  and  outer  0°  plies.  For  [90^ /On]*  composites,  tests  have  been  conducted  by 
Stm  et  al.  [46]  and  the  model  was  verified  in  the  previous  section. 

T300/976  Graphite/Epoxy  prepregs  were  selected  for  fabricating  both  flat  zmd 
cylindrical  test  panels,  which  were  cured  under  a  standard  cure  cycle  given  in  [32]. 
The  panels  were  then  sliced  into  specimens.  The  dimensions  and  the  ply  orienta¬ 
tions  of  the  specimens  are  listed  in  Table  3-2.  The  material  properties  for  Fiberite 
T300/976  are  shown  in  Table  3-3.  Each  specimen  was  X-rayed  to  examine  any  inter¬ 
nal  damage  resulting  from  manufacturing  or  cutting.  No  apparent  sign  of  damage 
was  found  in  the  X-radiographs  of  the  test  specimens. 

Each  specimen  was  clamped  firmly  at  both  ends  during  the  test.  The  fixtures 
fabricated  specifically  for  testing  of  flat  and  curved  specimens  are  shown  in  Fig¬ 
ure  3.12.  A  line- nosed  cylindrical  load  head  was  applied  downweird  at  the  centerline 
of  the  specimen.  The  configuration  of  the  test  setup  is  also  given  in  Figure  3.12. 
An  MTS  machine  was  utilized  for  all  the  tests.  The  load  head  speed  wais  0.0006 
in/min  based  on  displacement- controlled  loading.  A  strain  gage  was  mounted  at 
the  center  of  each  specimen.  The  strain  gage  reading  or  the  displacement  reading  as 
a  fimction  of  the  applied  load  was  recorded  during  the  entire  test  until  the  applied 
load  dropped  substantially,  as  a  result  of  internal  damage  in  the  specimens. 
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Figure  3.12  Schematics  of  the  testing  fixtures.  (Above);  The  testing  fixture  for 
flat  panels.  (Below);  The  testing  fixture  for  cylindriczd  panels. 
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Table  3-2 

Ply  orientations  and  geometry  of  test  specimens. 


Flat  Panels 

L  (in) 

W  (in) 

ff  (in) 

No.  of  Specimens 

[O6/9O3], 

1.5 

0.955 

0.098 

4 

[O4/9O4], 

1.5 

1.030 

0.0896 

4 

[O4/9O3/O2], 

1.5 

0.920 

0.098 

5 

[O3/9O2/O3/9O], 

1.5 

0.930 

0.098 

5 

Curved  Panels 

2a  (degree) 

W(in) 

H(in) 

No.  of  Specimens 

[O6/9O3], 

81 

0.970 

0.098 

5 

[08/902]j 

81 

0.930 

0.110 

5 

L  represents  the  span  length.  W  represents  the  plate  width.  H  represents  the 
plate  thickness.  2a  represents  the  arc  angle  for  curved  panels.  The  radius  for  the 
curved  panels  is  2.86  in. 
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Table  3>3 

Material  properties  for  Fiberite  T300/976. 


Material  Property 

Value 

Units 

Longitudinal  Young’s  modulus, 

17.6 

msi 

Transverse  Young’s  modulus,  Ey 

1.41 

msi 

Out-of-plane  Young’s  modulus,  Ez 

1.41 

msi 

Poisson’s  ratio, 

0.29 

Poisson’s  ratio,  u^z 

0.29 

Poisson’s  ratio,  Uyz 

0.40 

In-plane  shear  modulus,  Ggy 

0.81 

msi 

Out-of-plane  shear  modulus,  Gxz 

0.81 

msi 

Out-of-plane  shear  modulus,  Gyz 

0.50 

msi 

Longitudinal  thermal  expansion  coefficient,  0*1 

0.3 

■Jiil?. 

in‘F 

Transverse  thermal  expzuision  coefficient,  Oyy 

16.0 

tlin 

in^F 

Longitudinal  tensile  strength,  Xt 

220 

ksi 

Longitudinal  compressive  strength,  Xc 

231 

ksi 

Transverse  tensile  strength,  Yt 

6.46 

ksi 

Transverse  compressive  strength,  Yc 

36.7 

ksi 

Longitudinal  shear  strength  (cross  ply),  Sep 

15.5 

ksi 

Critical  streiin  energy  release  rate  for  Mode  I,  Gj^ 

0.50 

Ibf/in 

Critical  strain  energy  release  rate  for  Mode  II,  Gjic 

1.80 

Ibf/in 
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Flat  Panels 

The  test  data  of  flat  panels  with  [Oe/OOs],  and  [O4/9O4],  ply  orientations  are 
shown  in  Figure  3.13  and  Figure  3.14.  The  open  circles  shown  in  these  figures  are 
the  test  data.  As  can  be  seen  from  the  figures,  the  material  responses  between  the 
applied  load  and  the  measured  strains  were  quite  linear  until  the  maximum  load 
was  reached,  after  which  the  load  dropped  sharply.  Below  the  maximum  load,  no 
visual  damage  was  observed  during  the  test.  However,  significant  damage  including 
multiple  matrix  cracks  and  delaminations  was  found  immediately  after  the  load 
drop,  indicating  that  the  damage  occurrence  and  growth  were  not  progressive,  but 
unstable  and  catastrophic.  Apparently,  once  the  damage  initiated,  it  immediately 
propagated  extensively  inside  the  material  and  resulted  in  cataistrophic  failure.  A 
photograph  of  a  close  view  of  the  side  of  a  deunaged  specimen  neeur  the  center  region 
is  presented  in  Figure  3.15. 

Numerical  simulations  of  the  test  results  are  also  shown  in  Figure  3.13  and 
Figure  3.14,  represented  in  the  same  figures  by  the  solid  lines.  The  calculated 
load-strain  relationships  matched  with  the  data  very  well  up  to  the  final  failure 
load.  Initial  fmlure  due  to  matrix  cracking  in  the  central  90°  plies  was  predicted 
at  the  failure  load  in  both  cztses.  The  calculated  load  dropped  significantly  due 
to  an  extensive  delamination  which  occurred  at  the  interface  between  the  bottom 
ply  group  and  the  central  90°  plies,  resulting  from  the  matrix  cracks.  Following 
the  calculation,  once  the  matrix  cracks  initiated,  they  immediately  produced  zui 
interface  delamination  which  propagated  instantly  to  the  boundaries.  The  numerical 
simulations  of  the  deformed  configurations  of  a  [O4/9O4],  composite  specimen  as  a 
function  of  the  applied  load  are  presented  in  sequence  in  Figure  3.16.  It  is  noted 
that  additional  matrix  cracking  wets  also  predicted  at  other  locations  within  the 
central  90°  plies  during  the  imstable  delamination  propagation.  The  calculated 
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load  and  the  lower  interface  delamination  versus  the  displacement  of  the  load  head 
were  also  calculated  and  are  presented  in  Figure  3.17.  Clearly,  the  lower  interface 
delamination  grew  immediately  to  the  boundeiries  directly  after  the  critical  load  was 
reached. 

Cylindrical  Panels 

Figure  3.18  and  Figure  3.19  show  respectively  the  load-deflection  relationships 
of  [Oe/OOa],  and  [08/902]«  cylindrical  panels  subjected  to  concentrated  centred  line 
loading.  Symbols  in  both  flgures  represent  the  test  data,  and  the  solid  lines  are 
the  predictions  beised  on  the  model.  Each  symbol  corresponds  to  a  single  specimen. 
Both  types  of  specimens  resulted  in  the  same  type  of  damage  mode:  a  pair  of  shear 
cracks  occurring  in  the  90  degree  plies,  leading  to  an  upper  and  lower  interface 
delamination  along  the  0/90  interfaces.  A  photograph  of  a  side  view  of  a  [Og/OOa], 
specimen  after  testing  is  presented  in  Figure  3.15.  The  type  of  damage  was  basically 
similar  to  that  observed  in  [0/90],  flat  composites.  However,  the  overall  extent  of  the 
deleimination  along  the  bottom  0/90  interface  in  the  cylindrical  panels  was  smaller 
than  that  of  corresponding  flat  panels. 

The  predicted  curves  shown  in  Figure  3.18  and  Figure  3.19  agreed  with  the  test 
data  very  well  from  the  initial  loading  to  final  collapse.  As  predicted,  the  dzunage 
was  initiated  from  matrix  craicking  in  90  degree  plies  inside  the  laminates  and  the 
interface  delaminations  were  generated  from  the  matrix  cracking.  The  predicted 
locations  of  the  initial  matrix  cracks  and  the  extent  of  the  delamination  along  the 
bottom  0/90  interface  in  the  cylindrical  panels  as  compared  to  test  measurements 
are  presented  in  Table  3-4. 

A  typical  deformed  configuration  of  a  segment  of  a  [Oe/OOa],  cylindrical  compos¬ 
ite  corresponding  to  the  final  failure  load  is  presented  in  Figure  3.20.  A  significant 
amount  of  slip  deformation  can  be  seen  between  the  upper  and  lower  surfaces  of  the 
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BACK  STRAIN  (^in/in) 


Figure  3.13  Compzuison  between  the  measvured  and  the  calculated  load-strain  re¬ 
lationship  of  a  T300/976  (Oe/OOj],  composite  beam  subjected  to  a 
concentrated  line  load. 
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BACK  STRAIN  (Hin/in) 


Figure  3.14  Comparison  between  the  measured  and  the  cedculated  load-strain  re¬ 
lationship  of  a  T300/976  [O4/9O4],  composite  beam  subjected  to  a 
concentrated  line  load. 


Figure  3.15  Photographs  of  a  sideview  of  a  tested  [O4  /904]a  flat  panel  and  a  [Oe/QOsls 
curved  panel  near  the  loading  area. 
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Figure  3.16  Numerical  simulations  of  the  deformed  configurations  of  a  T300/976 
[O4/9O4],  composite  beam  subjected  to  a  transverse  concentrated  line 
load  at  various  loading  stages. 


Figure  3.17  The  csilculated  response  of  a  T300/976  [04/904],  composite  beam  sub¬ 
jected  to  a  transverse  concentrated  line  load.  (Above):  The  calculated 
applied  load  as  a  function  of  the  load  head  displacement.  (Below):  The 
extent  of  the  delamination  induced  by  the  matrix  crack  as  a  fimction 
of  the  load  head  displacement. 
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Table  3-4 

Matrix  crack  locations  and  delamination  extensions  for  curved  pEinels. 

{R  =  2.86  in,  2a  =  81°) 

Layups  ac/o  (model)  Q!c/a(data)  a^/a  (model)  a(f/a(data) 
[O6/9O3],  5.3  8.57  42  50.5 

[O4/9O4],  6  6.7  46.7  53.3 

All  the  numbers  in  the  table  are  in  terms  of  percent,  ttc  is  the  half  axe  angle 
corresponding  to  the  critical  matrix  crack  location,  qj  the  half  arc  angle  corre¬ 
sponding  to  the  final  delamination  extension  location,  and  a  the  half  arc  £mgle  for 
the  curved  panel. 
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delaxnination  at  the  lower  interface  delamination.  Again,  the  delamination  prop¬ 
agation  was  very  unstable.  The  delamination  immediately  propagated  away  from 
the  loading  area  and  led  to  the  collapse  of  the  structure.  It  was  calculated  that  the 
Mode  II  shear  fracture  dominated  the  deleimination  propagation  in  cylindrical  pan¬ 
els.  Accordingly,  it  seems  that  the  strain  energy  release  rate  Guc  is  more  criticed 
in  cylindrical  panels  than  in  flat  panels  for  governing  delamination  propagation. 

§3.5  Discussion 

Based  on  the  extensive  analytical  and  experimental  study,  matrix  cracking  is 
seen  to  be  the  initial  failure  mode  for  the  cross-ply  laminates  considered.  Delamina¬ 
tion  was  induced  by  these  matrix  cracks  along  the  interfaces  attached  to  the  matrix 
cracks.  For  90°  plies  on  the  bottom  surface,  a  matrix  crack  could  be  produced  due 
to  excessive  transverse  tensile  stresses  in  the  90®  plies.  These  types  of  cracks  axe 
referred  to  as  bending  cracks.  For  the  laminates  containing  90®  inner  plies,  ma¬ 
trix  failure  was  due  primarily  due  to  transverse  interlaminar  shear  stresses,  which 
reached  the  maximun  at  the  midplane  of  the  laminate.  This  type  of  inner  matrix 
cracking  is  referred  to  as  a  shear  crack.  The  location  of  the  shear  cracks  is  always  a 
distance  away  from  the  loading  point,  whereas  the  bending  crack  appears  directly 
beneath  the  loading  point. 

Accordingly,  depending  upon  the  sequence  of  the  ply  orientation  and  the  thick¬ 
ness  of  the  laminate,  either  bending  or  shear  cracks  could  occur,  earlier  than  the 
other.  For  thick  laminates,  interlaminar  shear  stresses  could  dominate  the  stress 
field,  leading  to  shear  cracking.  However,  for  thin  laminates,  bending  cracks  could 
appear  earlier  because  of  the  bending  deformation  of  the  plate.  The  growth  of  the 
delamination  induced  by  the  shear  cracks  is  apparently  very  unstable  and  catas¬ 
trophic,  whereas  the  bending  crack  induced  delamination  grows  stably  and  progres- 
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DISPLACEMENT  A  (in) 


Figure  3.18  Comparison  between  the  measured  and  the  calculated  load-strain  re¬ 
lationship  of  a  T300/976  [Oe/QOa],  curved  composite  beam  subjected 


to  a  concentrated  line  load. 


APPLIED  LOAD  PAV  (Ibs/in) 
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DISPLACEMENT  A  (in) 


Figure  3.19  Comparison  between  the  measured  and  the  calculated  load-strain  re¬ 
lationship  of  a  T300/976  [Og/OOj]*  curved  composite  beam  subjected 
to  a  concentrated  line  load. 


Chapter  3:  Damage  Induced  by  a  Cylindrical  Indenter 


55 


T300/976,  (0<i/903]s,  R  =  2.86  in,  2a  =  81° 


Figure  3.20  The  calculated  deformed  configuration  of  a  T300/976  [Oe/QOa],  curved 
composite  beam  corresponding  to  final  failure  load. 


FAILURE  LOAD  PAV  (Ibs/in) 
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STACKING  SEQUENCE 


Figure  3.21  The  failure  load  of  laminated  composite  be2uns  as  a  function  of  stack¬ 
ing  sequence.  Comparisons  between  the  measurements  amd  the  pre¬ 
dictions. 


FAILURE  LOAD  P/W  (Ibs/in) 
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Figure  3.22  The  predicted  failure  load  of  curved  laminated  beams  as  a  function  of 


curvature. 
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sively,  ajid  is  proportional  to  the  applied  load.  The  calculations  show  that  mixed 
mode  fracture  dominates  the  onset  of  the  shear  crack  induced  delamination  growth 
while  Mode  I  fracture  dominates  the  onset  of  the  bending  crack  induced  delami¬ 
nation.  Once  the  delamination  propagates,  Mode  I  fracture  continues  to  control 
the  growth  of  bending  crack-induced  delaminations,  but  Mode  II  fracture  becomes 
increasingly  important  for  the  growth  of  the  shear  crack  induced  delaminations. 

Apparently,  the  damage  mode  resulting  from  transverse  concentrated  line  loads 
strongly  depends  on  the  ply  orientation,  but  it  is  not  sensitive  to  the  change  of 
curvature.  However,  the  response  and  the  extent  of  damage  of  the  composites  are 
sensitive  to  both  ply  orientation  and  geometry  of  the  structures.  For  instance, 
Figure  3.21  shows  the  predicted  initial  failiire  load  corresponding  to  the  initiation 
of  matrix  cracking  for  different  stacking  sequences  of  flat  laminates.  Figure  3.22 
shows  the  effect  of  curvature  on  the  initial  and  final  failure  loads  of  cylindriceJ 
[Oe/OOa],  composite  panels.  Clearly,  ply  orientation  influences  the  initial  and  final 
failure  loads  of  the  structures  subjected  to  transversely  concentrated  loading. 

§3.6  Concluding  Remarks 

In  this  study,  a  2-D  zmalytical  model  has  been  developed  for  simulating  the 
response  of  laminated  composites  subjected  to  quasi-static  transverse  concentrated 
line  loads.  Experiments  were  also  performed  to  verify  the  einedysis.  Based  on 
the  study,  the  following  remarks  can  be  made  for  the  damage  growth  in  cross-ply 
composites  resulting  from  transverse  concentrated  line  loads: 

1.  Intraply  matrix  cracking  initiates  the  damage  in  laminated  composites. 

2.  Intraply  matrix  cracking  triggers  delamination. 

3.  Delamination  growth  induced  by  intraply  bending  cracks  is  stable  and  progres¬ 


sive. 
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4.  Delamination  growth  induced  by  intraply  shear  cracks  is  very  unstable  and 
catastrophic. 

5.  Mixed  mode  fracture  dominates  the  onset  of  shear  crack  induced  delamination 
while  Mode  I  fractxrre  dominates  the  onset  of  bending  crack  induced  delamina¬ 
tion. 

6.  The  growth  of  the  delamination  induced  by  a  bending  crack  is  governed  by 
Mode  I  fractiire. 

7.  The  growth  of  the  delamination  induced  by  sheair  cracks  depends  strongly  on 
the  Gjic  value. 

8.  Both  stacking  sequence  and  curvature  affect  the  response  and  the  extent  of  the 
d2Lmage  in  composites. 


Chapter  4 


Damage  Induced  by  a  Spherical  Indenter 


§4.1  Statement  of  Problem 

The  objective  of  this  study  was  to  fundamentally  understand  delamination 
propagation  behavior  and  the  effect  of  matrix  cracking  on  its  growth  in  cross-ply 
composites  resulting  from  a  transverse  concentrated  load.  Therefore,  only  [0m/90n]* 
composite  plates  were  considered.  The  rectangular  plates  could  be  simply  supported 
or  clamped  on  any  one  of  the  edges  and  subjected  to  a  transverse  quasi-static 
concentrated  load  by  a  sphericzd  indenter. 

In  order  to  achieve  the  objective,  the  plates  were  assumed  to  contain  a  pre¬ 
existing  small  delamination  located  at  the  centred  loading  area  on  the  second  90/0 
interface  directly  beneath  the  indenter.  Two  types  of  matrix  cracking  in  conjimction 
with  the  delamination  were  considered  in  the  study:  a  bending  crack  emd  a  pair  of 
shear  cracks.  Accordingly,  for  a  given  loading  and  boimdary  condition,  four  types 
of  interned  deunage  modes  were  studied:  1)  deleunination  with  no  matrix  cracks;  2) 
deleunination  induced  by  a  bending  crack;  3)  deleunination  induced  by  sheeu^  creicks; 
and  4)  delamination  induced  by  both  a  bending  crack  aind  shear  cracks,  as  shown 
in  Figure  4.1. 
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Figure  4.1  Four  types  of  damage  modes  considered  in  this  einalysis:  1)  delami¬ 
nation  with  no  matrix  cracks;  2)  deleunination  induced  by  a  bending 
crack;  3)  delamination  induced  by  shear  cracks;  and  4)  delamination 
induced  by  both  a  bending  crack  emd  shear  cracks. 
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The  analytical  model  proposed  for  describing  the  behavior  of  laminated  com¬ 
posite  panels  consists  of  three  portions:  a  stress  analysis,  a  contact  analysis,  and  a 
failure  analysis.  The  stress  ancdysis  is  based  on  the  large  deformation  theory  sim¬ 
ilar  to  the  one  presented  in  Chapter  3,  except  that  the  displacement  field  is  now 
three-dimensional.  Therefore,  the  theory  will  not  be  elaborated  in  this  chapter. 

4.2.1  Contact  Analysis 

Two  types  of  contact  axe  involved:  the  rigid-elastic  contact  between  the  inden¬ 
ter  and  the  plate,  and  the  elastic-elastic  contact  between  the  cracked/delaminated 
interfaces,  as  shown  in  Figure  4.2.  The  local  indentation  resulting  from  a  spheri¬ 
cal  indenter  is  very  complicated  and  three-dimensional.  The  local  indentation  could 
significantly  affect  the  delamination  growth  and  its  interface  deformation,  especially 
when  the  delamination  is  small.  The  actual  contact  between  the  sphericeil  indenter 
and  the  laminate  had  to  be  simulated.  For  simplicity,  only  the  elzistic-elastic  contact 
is  presented  in  this  section.  The  rigid-elastic  contact  can  be  treated  as  a  special 
case  of  the  elastic-eleistic  contact.  The  contributions  of  more  thEin  one  contact  pair 
can  be  easily  handled  as  the  formulation  is  genereJ  for  contact  between  multiple 
bodies. 

Several  methods  exist  to  implement  the  contact  constraints.  The  most  popular 
choices  are  the  Lagrange  multiplier  method  emd  the  penalty  method.  The  former 
has  the  advantage  of  enforcing  the  exact  constraints,  but  induces  additional  param¬ 
eters  which  cause  an  increase  of  the  finite  element  matrix  bandwidth  and  therefore 
substantially  enlarge  the  overall  size  of  the  equation  system  to  be  solved  [56].  This 
method  was  adopted  for  the  study  of  damage  induced  by  a  cylindrical  indenter  for 
analyzing  the  interaction  between  matrix  cracking  and  delamination.  However,  the 
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Figure  4.2  Contact  problem  associated  with  a  delaminated  leoninate  subjected  to 


an  indenter. 
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application  of  this  method  to  three  dimensions  becomes  extremely  costly.  Although 
the  penalty  method  has  the  advantage  of  not  requiring  additional  equations,  it  is 
very  sensitive  to  the  choice  of  a  penalty  parameter,  which  can  possibly  lead  to 
ill-conditioning  when  the  value  of  penalty  parameter  increases  [56,  62,  66]. 

Therefore,  in  this  study,  an  augmented  Lagrangian  method  was  adopted.  It 
was  originally  proposed  by  Hestenes  [57]  and  Powell  [58]  in  studying  mathemati¬ 
cal  programming  problems  subject  to  equality  constraints.  It  has  been  shown  to 
provide  important  advantages  over  the  more  traditional  Lagr2mge  multiplier  and 
penalty  methods.  The  applicability  of  this  method  to  mathematical  progreunming 
problems  subject  to  inequality  constraints  is  also  well  documented  by  Rockafellar 
[59]  and  others.  Detailed  mathematical  discussions  can  be  found  in  [56].  More  re¬ 
cently,  this  technique  has  been  successfully  applied  to  various  constrained  mechanics 
problems  such  as  incompressible  finite  deformation  elasticity  [60],  and  contact  prob¬ 
lems  in  two  dimensions  [61,  66].  The  augmented  Lagrangian  techniques  have  been 
known  to  provide  almost  exact  enforcement  of  constraints  while  using  finite  penalty 
parameters,  avoiding  the  ill-conditioning. 

Using  the  Lagrange  multiplier  technique,  the  total  potential  energy  with  the 
contact  constraint  conditions  can  be  described  as 


n  =  n-h 


(4.1) 


The  variational  equations  of  the  laminate  subjected  to  a  spherical  indenter, 
corresponding  to  the  Lagrange  multiplier  formulation  (similar  to  Eqs.  (3.8)),  can 
be  described  as 


^n-t- 


L 

I 


XN^gdT  =  0, 

SXpjg  dT  =  0 


1 


(4.2) 


where  A  at  is  unknown,  which  needs  to  be  solved  from  the  above  equation.  The 
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introduction  of  the  Lagrange  mtiltiplier  \n  could  considerably  enlarge  the  system 
of  equations. 

However,  the  penalty  treatment  is  achieved  by  the  replacement  of  A  a/  in  Eqs.  (4.2) 
by  esg,  where  is  defined  as  the  penalty  parameter.  It  can  be  observed  that  as 
€n  oo,  g  —*  0  and  A^v  is  bounded.  However,  as  es  —*  oo,  the  gap  g  must  be  zero, 
which  satisfies  the  constraint  condition. 

The  variational  equation  for  the  penalty  method  can  now  be  written  as 

<5n  +  y*  €NgSgdr  =  0  (4.3) 

It  is  noted  that  Eq.  (4.3)  only  involves  the  displacement  variables,  emd  no  addi¬ 
tional  contact  force  constraint  is  needed.  As  a  result,  Eq.  (4.3)  does  not  introduce 
additional  equations  and  is  extremely  attractive  for  finite  element  implementations. 
However,  the  constraint  condition  g  >  0  is  only  satisfied  when  approaches  in¬ 
finity,  which  unfortunately  leads  to  ill-conditioning.  In  order  to  solve  Elq.  (4.3), 
eN  is  chosen  as  a  large  finite  constant  in  practice  without  inducing  ill-conditioning. 
Therefore,  some  overlapping  between  contact  surfaces  may  occur.  In  the  present  de¬ 
lamination  contact  problem,  the  overlapping  between  the  contact  surfaces  needed 
to  be  minimized.  These  considerations  led  to  the  development  of  the  method  of 
augmented  Lagrangians. 

In  the  method  of  augmented  Lagrangians,  the  Xn  is  initially  chosen  to  be 
an  arbitrary  constant.  If  this  constamt  is  not  the  correct  Lagrange  multiplier,  the 
contact  constraint  is  not  satisfied  and  minimization  of  the  total  potential  energy  of 
Eq.  (4.1)  does  not  lead  to  the  equation  of  equilibrium.  Therefore,  from  a  penalty 
function  viewpoint,  the  total  potential  energy  represented  by  Eq.  (4.1)  needs  to  be 
further  penalized  by  the  following  modified  functional: 
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n  =  n  +  dr  +  ^  dr  (4.4) 

where  <  0  denotes  the  fixed  estimate  of  the  correct  Xfj.  The  superscript  k 
indicates  that  the  search  for  the  correct  A is  an  iterative  process.  The  updated 
formula  is 

A<„*+"  =  +  e^g)  (4.5) 

Therefore,  the  variational  equation  corresponding  to  Eq.  (4.4)  can  be  derived 
as 


OT  +  (A^;^^  +  eN9)  SgdT  =  0  (4.6) 

It  is  noted  that  the  term  (A^^  +  cnq)  plays  the  role  of  the  Lagrange  multiplier 
A/v-  If  A;v  is  the  correct  multiplier,  then  g  =  0  on  T.  Thus,  in  the  case  where 
the  multipliers  are  correct,  Eq.  (4.6)  achieves  exactly  the  same  form  eis  the  first  of 
Eqs.  (4.2),  making  it  an  exact  peneJization. 

It  is  important  to  notice  that  Eq.  (4.6)  is  a  nonlinear  equation  due  to  the 
contact  conditions  and  geometric  nonlinearity.  In  general,  then,  it  will  be  necesseury 
to  solve  Eq.  (4.6)  in  an  iterative  manner.  In  practice,  is  chosen  to  be  eis  large  as 
practically  possible  without  inducing  ill-conditioning.  The  advantage  of  the  current 
treatment  over  the  penalty  method  is  that  satisfaction  of  the  constraints  can  be 
improved  even  if  is  of  relatively  modest  value  through  repeated  application  of 
the  augmentation  procedure.  Since  these  augmentations  only  chzinge  which 
is  fixed  with  regard  to  solution  of  Eq.  (4.6),  the  ill-conditioning  problem  usually 
associated  with  the  penalty  method  is  mediated  or  eliminated. 
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The  initiation  of  delamination  growth  (onset  of  delamination  growth)  was  pre¬ 
dicted  based  on  lineax  elastic  fracture  mechanics.  The  well-known  virtual  crack 
closure  technique  [52]  served  as  the  basis  for  the  strain-energy  release  rate  calcu¬ 
lation.  This  procedure  determines  Mode  I,  Mode  II,  and  Mode  III  strain  energy 
release  rates  from  the  energy  required  to  close  the  delamination  over  a  small  area. 
The  strain  energy  release  rates  Gj,  Gu,  and  Gm  for  Mode  I,  Mode  II,  Jind  Mode 
III  fracture  respectively,  can  be  expressed  as  [75] 


Gj  =  lim  i 

AA->o  \ 

to 

> 

fAA 

/  [u;['(A-|- AA)-ii“(A-|- AA)]  <T„„dA 

0 

1  (4.7.a) 

G//  =  lim  \ 
AA-^o  I 

f  1 

[2AA. 

fAA  ] 

j  [u+(A-|-AA)-u7(A-|-AA)]a,„dAj 

1  (4.7.!.) 

Gill  =  lim 

AA~*o  1 

f  1 

[2AA, 

fAA  ] 

j  [u^(A -f  AA)  —  Uj  (A -1- AA)]  (Ttn  dA 
Jo  j 

[  (4.7.C) 

where  AA  is  the  crack  extension  eurea;  n,  s,  /  form  a  local  coordinate  system  which 
axe  out-of-plane  normal,  in-plane  normal,  and  in-pleine  tiingent  respectively;  cr„„, 
(Tan,  and  (Ttn  are  the  stresses  at  the  crack  front  associated  with  a  crack  size  of  A;  and 
U+,  u^,  u/"  and  u“,  uj,  uj"  are  the  displacements  of  the  upper  and  lower  svufaces 
of  the  crack  eissociated  with  a  crack  size  of  A  -I-  AA  respectively. 

Therefore,  at  any  point  on  the  delamination  front,  strain  energy  release  rates 
can  be  calculated.  Since  the  delamination  growth  can  be  attributed  to  a  mixed  mode 
fractiire,  the  criterion  for  determining  the  initiation  of  the  delamination  growth  was 
selected  eis  [75-78] 


for  any  point  on  the  delamination  front.  Here,  ^///c 

critical  strain  energy  release  rates  corresponding  to  Mode  I,  Mode  II,  and  Mode  III 


Chapter  4:  Damage  Induced  by  a  Spherical  Indenter 


68 


fracture,  respectively.  It  was  assumed  that  G/c  G//,.,  £ind  G///^  did  not  change 
with  delamination  size.  Based  on  the  two-dimensional  study  in  this  thesis  and  [42], 
a  =  1,  /3  =  1,  and  7  =  1  were  selected  for  this  study,  because  they  were  found 
to  provide  the  best  fit  to  the  experiments  in  two  dimensions.  It  was  assumed  that 
GjUc  —  Gjic  [67],  because  the  value  of  Gjjj  has  not  been  reported  in  the  literature. 
Accordingly,  delamination  would  start  to  propagate  when  >  1  at  any  point  on 
the  delamination  front. 

§4.3  Finite  Element  Implementation  of  the  Formulation 

The  laminate  was  discretized  by  8-node  isopau-ametric  brick  elements  into  n*/ 
element  domains  with  n„p  nodal  points. 

The  displacement  field  within  an  element  was  assumed  to  be  in  the  form 

^en 

U  =  '^NaUa  (4.9) 

a=l 

where  Na  is  the  displacement  shape  (interpolation)  function  associated  with  the 
local  nodal  number  a,  and  Ua  is  the  nodal  displacement  vector  at  node  number  A. 
The  displacement  shape  functions  Na  in  the  local  coordinates  p,  g,  and  r  have  the 
following  expression: 

Naip,q,r)  =  i  (1  -f-ppa)(l  -|-gg„)  (1  -|-rra)  (4.10) 

where  a  =  1  —  8  and  there  is  no  summation  on  a. 

Both  the  indenter  and  the  laminate  were  discretized  by  finite  numbers  of  ele¬ 
ments,  and  the  master-slave  definition  is  the  same  as  given  in  the  2-D  study.  The 
indenter  was  treated  as  the  master  surface  while  the  top  of  the  laminate  was  treated 
as  the  slave  surface.  For  delamination  surfaces,  both  upper  zind  lower  surfaces  may 
be  treated  as  slave  and  master  surfaces.  It  was  again  assumed  that  the  contact  pres- 
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sure  on  each  master  element  was  constant.  Therefore,  Eq.  (3.21)  and  its  associated 
arguments  still  applied. 


4.3.1  Variational  Formulation  of  the  Discrete  Problem 

The  finite  element  formulation  employs  the  modified  total  potential  energy 
(Eq.  (4.4))  in  its  discrete  form.  Thefore,  the  following  fimctional  is  obtained 

n(u)  =  n(u)  +  (4.11) 

where  the  components  of  vector  A  are  assumed  to  be  a  fixed  estimate  of  the  correct 
Lagrange  multiplier  vector  and  need  to  be  updated  following  the  formula  of  Eq.  (4.5) 
in  a  vector  form.  The  focus  here  is  on  the  contribution  of  the  contact  constraints  due 
to  cracked/delaminated  interfaces.  The  discrete  form  of  the  virtual  work  principle 
can  be  derived  from  Eq.  (4.4)  and  has  the  following  form: 

(^+t^|6)iu  =  0  (4.12) 

imagined  as  the  approximate  contact  force  vector.  The  t  will  be  the  correct  contact 
force  vector  if  the  gap  vector  g  is  essentially  zero.  It  is  reiterated  that  A  is  a  constant 
vector  in  the  solution  of  Eq.  (4.12).  Therefore,  the  nonlinear  equilibrium  equations 
for  the  discrete  system  can  be  obtained  as 


du  du 


(4.13) 


In  the  framework  of  the  finite  element  method,  the  term  can  be  replaced  by 
Pint  —  Rexti  defined  in  the  same  way  as  in  the  2-D  study.  The  term  can  be 

replaced  by  Fc,  which  is  a  contact  force  vector.  Therefore  the  following  equation  is 
obtained: 


Pint  +  Fc  —  R«xt 


(4.14) 
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The  equation  represents  a  weak  form  of  the  equilibrium  equation  in  the  framework 
of  Lagrangian  augmentation  technique.  It  is  shown  that  t  represents  the  vector  of 
contact  forces,  which  are  distributed  according  to  the  matrix  to  the  associated 
nodes  of  the  contact  element. 

In  general,  due  to  the  finite  deformation  and  the  contact  constraints,  the  sys¬ 
tem  of  equations  is  nonlinear  and  therefore  needs  to  be  solved  iteratively  by  the 
Newton- Raphson  method  [70-74].  By  applying  the  Newton- Raphson  technique, 
the  following  linearized  equation  is  obtmned: 

a 

^[Pint  +  Fcjdu  =  Rext  —  (Pint  +  Fc)  (4-15) 

Recalling  that  Pint  and  Fc  are  the  function  of  u  only,  the  above  equation  can 
be  written  as 


l|;(Pint)  +  =  R.xt  -  (Pint  +  F.) 


(4.16) 


Or 


(K  Kc)du  =  iUxt  -  (Pint  +  Fe)  (4.17) 


where  K  (=  •;^(Pint))  is  the  cletssical  tangent  stiffness  matrix  which  can  be  obtained 
without  considering  the  contact  constraint  conditions.  Fc  is  the  contact  residueil 
forte  vector  and  defined  as 


Fc 


(4.18) 


and  Kc  is  the  contribution  of  the  contact  conditions  to  the  tangent  stiffness  matrix, 
which  guarantees  the  best  approximation  to  the  actued  response.  It  is  defined  as 
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For  convenience,  the  above  matrices  will  be  derived  on  the  basis  of  an  element 
coming  in  contact  with  one  node  on  the  other  body,  forming  one  contact  element. 
Such  an  element  provides  one  gap  g  and  one  estimated  contact  force  magnitude  t, 
and  has  N  nodal  points.  Subsequently,  capital  letter  indices  {N  and  M)  are  used 
here  to  identify  quantities  which  are  related  to  the  element  nodal  points.  Moreover, 
for  repeated  indices  the  summation  convention  of  tensor  calculus  holds.  The  contact 
residual  force  of  an  element  is 

t  (4.20) 

where  t  =  X  + eg  and  A  is  a  constant  Lagrange  multiplier  in  some  step.  The  stiffness 
contribution  due  to  contact  has  the  following  form 


,  dU  ^  fx  1  ^  r/x  \  ^g  1 


(4.21) 


ke  = 


(4.22) 


_  dg  dg  d  dg 

'  ”  auN  auM  ^  auN'^auM^ 


(4.23) 


The  corresponding  form  of  Eq.  (4.16)  for  one  contewrt  element  is 


(•  ^^int  I  ^9  ^9  ,  a  .  dg  \x  1  N  _  pM  _  /^pM  I 

^  au^  au^auA^  -*^ext  (Pmt  +  Fc ) 


(4.24) 


The  three  terms  in  parentheses  form  the  tangent  stiffness  matrix  of  the  element. 
The  first  term  is  the  classical  tangent  stiffness  matrix  for  undamaged  composites. 
The  derivation  of  the  matrix  is  available  in  the  literature  related  to  finite  element 
methods  [70-74]. 
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The  derivation  of  the  contact  element  stiffness  matrix  kc  is  important,  as  it  ren¬ 
ders  the  quadratic  convergence  of  Newton’s  method.  The  derivation  of  the  contact 
stiffness  matrix  kc  and  the  residual  force  fi;  are  described  in  Appendix  D. 

4.3.2  Mesh  Generation 

It  has  been  observed  experimentally  that  delaminations  in  cross-ply  composites 
due  to  low- velocity  impact  or  quasi-static  indentation  loading  form  ellipsis-like  or 
peanut-like  shapes.  Therefore,  only  delaminations  with  circular  or  elliptical  shapes 
were  considered  in  the  analysis,  as  shown  in  Figure  4.3,  Figure  4.4,  Figure  4.5,  and 
Figure  4.6. 

This  section  describes  the  technique  used  for  the  finite  element  mesh  generation 
for  the  proposed  problems.  The  procedure  used  for  generating  the  meshes  wets 
similar  to  the  one  given  by  Whitcomb  [68j.  The  meshes  by  Whitcomb  were  used  only 
for  modeling  a  delamination  in  a  square  laminate  without  including  matrix  cracking. 
The  present  meshes  work  for  a  rectanguleir  laminate  conteuning  a  deleimination  in 
conjunction  with  matrix  cracks.  The  procedure  of  constructing  a  typical  mesh  is 
described  as  follows. 

It  was  assumed  that  the  length  of  the  plate  is  greater  than  or  equal  to  its  width. 
First,  a  two-dimensionzd  mesh  was  formed  inY  —  Z  plane,  where  the  delamination 
is  represented  by  a  line  of  size  6.  Here,  b  is  the  size  of  the  short  ajcis  of  the  elliptical 
delamination.  Secondly,  this  two-dimensional  mesh  was  swept  through  the  space 
by  a  90°  arc  around  the  Z  axis  to  generate  a  1/4  cylindrical  3-D  mesh.  Here,  a 
circular  delamination  of  radius  b  was  obtained  in  this  cylindriced  mesh.  The  radial 
lines  were  normal  to  the  circumferential  lines.  Third,  the  outer  part  (it  cam  be 
more  than  one  surface)  was  then  transformed  to  obtain  a  square  boundary.  The 
fourth  step  was  that  the  inner  circle  was  replaced  by  a  block  consisting  of  brick 
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Crack  Length 


Figure  4.6  The  coordinate  system  used  in  the  finite  element  analysis  for  a  lam¬ 
inate  containing  a  pair  of  shear  cracks  and  one  delamination  having 
two  semi-circular  or  semi-elliptical  shapes  offset  by  the  internal  cracks. 
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elements  so  that  only  brick  elements  needed  to  be  used.  The  fifth  step  was  to  apply 
a  modified  elliptical  transformation  to  obtain  an  ellipticed  delamination  front,  while 
keeping  good  spacing  for  elements  near  the  center.  Based  on  the  standaird  elliptical 
transformation,  a  scale  factor  was  introduced  for  “stretching”  the  mesh,  for  example, 
in  longitudinal  or  width  (AT)  direction: 

=  (4.25) 

where 

./Y2  I  y2  _ _ 

F  =  - - - ,  /or  r  >  >/A-2  +  (4.26) 

T 

and 

F  =  1,  /or  r  <  (4.27) 

The  mesh  can  be  similarly  “stretched”  in  Y  direction  if  6  >  a.  It  is  easy  to  show 
that  when  F  =  1,  the  transformation  is  the  standard  elliptical  (conformal)  mapping. 
In  the  conformal  mapping,  any  two  lines  which  were  orthogon2il  would  be  kept 
orthogonal  after  the  mapping.  When  F  <  1  the  mapping  is  a  modified  elliptical 
mapping,  which  can  produce  a  more  evenly  distributed  mesh  for  the  mesh  within  a 
circle  of  radius  of  r.  By  choosing  the  parameter  r  a  little  less  than  the  radius  of  the 
delamination  front  in  the  cylindrical  mesh,  the  orthogonality  was  well  maintained 
in  the  neighborhood  of  the  delamination  front  in  the  ellipticeJ  trzinsformation.  By 
maintaining  orthogonality,  the  calculation  of  fracture  mechanics  quantities  such  as 
energy  release  rates  can  be  greatly  simplified. 

For  modeling  the  surface  crack,  it  was  assumed  that  the  nodes  on  the  crack 
surface  were  free  to  move  and  had  independent  degrees  of  freedom.  Due  to  symme¬ 
try,  only  1/4  laminate  was  modeled  and  therefore  only  one  sheair  crack  needed  to 
be  modeled.  In  meshing  the  plate  with  a  shear  crack,  the  laminate  was  divided  into 
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two  blocks,  the  delaminated  block  and  the  solid  block,  separated  by  the  shezir  crack 
surfaces.  The  mesh  of  the  delaminated  block  away  from  the  indenter  was  formed 
first,  by  the  procedure  described  above.  Then,  the  mesh  of  the  solid  block  which 
was  below  or  near  the  indenter  was  generated  and  attached  to  the  mesh  of  the  first 
block.  The  nodes  on  both  sides  of  the  shear  crack  were  assumed  to  move  freely.  A 
typical  mesh  used  forthe  analysis  proposed  is  shown  in  Figure  4.7. 

4.3.3  Contact  Node  Search 

The  general  approach  to  contact  in  the  literature  is  the  slave-master  method. 
That  is,  one  side  of  the  contact  surface  is  designated  as  the  master  surface  while  the 
other  is  treated  as  the  slave  surface.  In  the  finite  element  method,  this  practice  is 
referred  to  as  one-pass  and  obviously  lacks  symmetry.  In  the  current  delamination 
problem  where  comers  and  edges  exist,  and  where  complicated  curvature  also  exists 
during  deformation,  the  lack  of  symmetry  can  produce  slow  convergence  or  diver¬ 
gence  and  incorrect  solutions.  The  remedy  may  be  easily  achieved  by  periodiczdly 
reversing  the  designation  of  the  master  and  slave  surfaces.  A  more  sophisticated 
approach  is  adopted  here  which  makes  two  passes  at  every  time  step  with  each 
surface  being  the  master  surface  on  one  pass  [64].  Taking  two  passes  each  load  step 
doubles  the  cost  of  the  algorithm.  However  for  many  practical  problems,  a  sym¬ 
metric  zdgorithm  is  much  more  robust  than  a  conventional  unsymmetric  algorithm, 
making  it  well  worth  the  additional  cost. 

In  the  two-dimensional  study  of  this  thesis,  a  structured  finite  element  mesh 
was  used,  which  required  less  bookkeeping  in  terms  of  the  contact  search.  However, 
in  the  current  three-dimensional  study,  an  unstructured  finite  element  mesh  was 
used.  It  vfas  very  challenging  to  find  an  economical  search  algorithm. 

In  a  finite  element  setting,  given  the  slave  node,  the  problem  remaining  is  to 


Figure  4.7  Description  of  a  typical  finite  element  mesh  used  in  the  analysis. 
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seeirch  for  the  closest  projection  onto  the  master  surface.  There  are  three  steps  to 
take  for  the  search.  The  first  step  is  to  determine  which  master  node  is  closest 
to  the  slave  node.  The  second  step  is  to  determine  which  of  the  master  segments, 
having  the  master  node  as  one  of  their  vertices,  is  closest  to  the  slave  node.  The 
third  step  is  to  determine  the  closest  point  of  the  master  segment  to  the  slave  node 
in  terms  of  ^  and  rj  in  the  isoparametric  biunit  domain.  For  each  slave  node,  one 
contact  point  can  be  located.  Therefore,  the  gap  g  can  be  determined.  Due  to  the 
complexity  of  the  contact  search  process,  the  details  of  these  steps  are  deferred  to 
Appendix  C. 
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The  delamination  growth  model  was  implemented  into  a  finite  element  analysis. 
The  modified  crack  closure  technique  [52]  served  as  the  basis  for  the  strain  energy 
relezise  rate  calculation.  The  physical  interpretation  of  the  crack  closure  technique  is 
that  the  amount  of  work  required  to  close  the  crack  by  a  short  distance  determines 
the  energy  release  rates  Gi,  G//,  and  Giij.  The  crack  closure  energy  involves 
products  of  forces  and  displacements  of  elements  along  the  delamination  front. 

The  energy  release  rate  calculation  is  illustrated  for  an  8-node  brick  element 
in  Figure  4.8  which  shows  a  schematic  of  the  delamination  front  region.  The  nodes 
of  interest  for  the  energy  release  rate  calculations  are  indicated  by  the  solid  circles. 
In  the  finite  element  mesh,  it  is  advantegeous  to  use  a  more  regular  mesh  near  the 
delamination  front  so  that  a  single  solution  method  can  be  used  [52]. 

In  calculating  the  crack  closure  work,  it  was  necessary  to  calculate  the  forces 
transmitted  at  the  crack  tip  as  well  as  the  crack  front  relative  displacements.  The 
forces  were  evzJuated  by  summing  up  all  the  nodal  forces  for  all  the  elements  which 
were  connected  to  the  nodes  on  the  front  and  whose  centroids  lay  above  the  delam¬ 
ination  plzuie. 

Energy  release  rate  is  a  measure  of  the  energy  per  unit  area.  Hence,  the  work 
must  be  normalized  by  the  appropriate  areas.  Unfortunately,  there  is  no  simple 
or  exact  way  to  determine  the  appropriate  areas.  Certmn  approximations  must  be 
m£wle.  It  was  assumed  that,  for  the  linear  elements,  the  work  associated  with  the 
elements  attached  to  line  BE  (see  Figure  4.8)  was  proportional  to  their  respective 
areas.  For  radial  line  BE,  the  area  associated  with  it  was  estimated  as  follows: 

AA  =  ^{Ar+A2)  (4.28) 

where  the  physical  interpretation  was  that  the  weights  for  the  work  produced  by 
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each  element  were  proportional  to  its  area. 

As  the  delamination  front  in  general  was  not  parallel  to  one  of  the  global 
coordinate  axes,  a  local  orthogonal  coordinate  system  was  defined  in  which  the 
forces  and  relative  displacements  could  be  evaluated  in  terms  of  the  local  coordinates 
see  Figure  4.8).  The  local  coordinate  system  had  one  axis  tangent  to  the 
delamination  front,  one  axis  normal  to  the  delamination  front,  and  one  axis  normal 
to  the  delamination  plane.  The  transformed  nodal  force  vector  in  the  local  system 
can  be  evaluated  as 

F'i  =  TyFj  (4.29) 

where  T  is  the  transformation  matrix  of  rank  3x3  from  the  global  coordinate  system 
to  local  coordinate  system.  And  the  components  of  T  are  directional  cosines  of  6 
(shown  in  Figure  4.8). 

Therefore,  the  energy  release  rates  for  Modes  I,  II,  and  III  can  be  calculated  in 
one  step: 

(4.30.i) 

Once  the  energy  release  rates  were  ctdculated,  the  delamination  growth  criterion 
was  applied  for  each  point  on  the  delamination  front.  In  order  to  predict  the 
delamination  growth,  the  three-dimensional  crack  growth  criterion  (Eq.  (4.8))  was 
used.  U  Ed  >  1,  then  delamination  growth  was  predicted. 
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The  proposed  3-D  FEM  model  was  verified  by  comparing  calculations  with 
existing  analytical  and  test  results. 

A  DCB  under  a  Line  Load 

It  was  shown  by  Davidson  ei  al.  [81-83]  that  a  DCB  I  specimen  with  a  straight 
delamination  front  had  the  highest  value  of  strain  energy  release  rate  Gi  at  the 
center  and  lowest  value  at  its  edges.  The  solutions  were  obtained  in  three  dimen¬ 
sions,  using  the  commercial  finite  element  code  ABAQUS.  The  material  used  was 
Hexcel  T300/F155  graphite  epoxy.  Its  properties  are  shown  in  Table  4-1.  Similar 
numerical  Ccdculations  were  performed  using  the  current  code,  and  the  results  were 
compared  with  the  ones  given  by  Davidson  et  al.  [81]. 


Table  4-1 

Material  properties  for  Fiberite  T300/F155. 


Material  Property 

Value 

Units 

Longitudinal  Yoimg’s  modulus.  Ex 

17.9 

msi 

Transverse  Young’s  modulus,  Ey 

1.105 

msi 

Out-of-pleine  Young’s  modulus,  Ez 

1.105 

msi 

Poisson’s  ratio,  Uxy 

0.34 

Poisson’s  ratio,  i/j. 

0.34 

Poisson’s  ratio,  Uyz 

0.34 

In-plane  shear  modulus,  Gxy 

0.5 

msi 

Out-of-plane  shear  modulus,  Gn 

0.5 

msi 

Out-of-plauie  shear  modulus,  Gyz 

0.5 

msi 
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For  DCB  in  Mode  I  (see  Figure  4.9),  the  comparison  between  the  solutions 
from  [83]  zind  the  current  code  agreed  quite  well.  The  G/  in  the  middle  portion  of 
the  laminate  for  [0)24  was  fairly  flat,  but  it  was  very  sensitive  to  the  ply  orientation. 
For  [0)24,  the  Gj  value  varies  signiflcantly  along  the  width  of  the  specimen.  It  is 
believed  that  this  phenomenon  occurs  because,  when  the  beam  is  in  bending,  the 
upper  and  lower  arms  of  the  beam  deform  anticlastically  and  tends  to  close  the 
crack  near  the  free  edge.  As  a  result,  the  crack  surfaces  neEir  the  free  edge  may 
be  in  contact  due  to  the  anticlastic  curvature  effect  for  emgle  ply  leuninates.  For 
[45/  —  45/(— 45/45)2]2s,  it  was  foimd  that  the  finite  element  nodes  on  upper  and 
lower  arms  of  the  DCB  specimen  penetrated  into  each  other  without  imposing  the 
contact  model,  as  shown  in  Figure  4.9.  By  using  the  contact  model,  the  overlapping 
between  the  edges  of  the  upper  and  lower  arms  was  prevented,  resulting  in  the  zero 
Gi  value  at  the  edge  nodes.  The  observation  that  a  smadl  contact  zirea  at  the  free 
edge  of  the  crack  front  actually  agreed  well  with  the  experimental  findings  [81,  82]. 

Therefore,  the  G/  was  the  highest  in  the  middle  of  the  laminates,  and  the 
middle  portion  would  grow  earlier  when  the  DCB  specimens  were  subjected  to  a 
line  load.  Actually,  it  was  observed  by  Davidson  ei  al.  [81]  that  in  DCB  Mode  I 
testing,  the  delamination  front  was  somewhat  thumbnail  shaped. 

Numerical  simulations  were  also  performed  on  the  DCB  specimens  due  to 
Mode  II  loading.  The  solutions  show  that  the  G//  values  were  fairly  constant 
for  the  [0]24  laminate,  but  varies  very  drastically  along  the  width  of  the  specimen 
for  angle  ply  laminate.  The  G//  distribution  wais  very  much  different  from  that 
of  the  Mode  I  for  the  same  ply  orientation.  Surprisingly,  it  was  found  from  the 
calculation  that  only  a  narrow  strip  of  elements  were  in  contact  along  the  loading 
edge,  while  the  rest  of  the  crack  surfaces  was  open. 
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Figure  4.9  Distribution  of  Gj  for  angle-ply  DCB  specimens  with  strzught  delam 


ination  fronts  in  Mode  I.  (W  =  1  in,  c  ^  1  in.) 
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Figure  4.10  Distribution  of  Gjj  for  angle-ply  DCB  specimens  with  straight  delam¬ 
ination  fronts  in  Mode  II.  (IF  =  1  in,  c  =  1  in.) 
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A  closed  form  solution  (CFS)  for  a  circular  isotropic  plate  due  to  a  central  point 
load  is  given  in  [80].  The  solution  is  exact  for  linear  deflections  and  approximate 
for  large  deflections.  The  transverse  deflections  at  the  center  of  the  plate  were 
determined  by  the  cxurent  code  as  a  fxmction  of  the  applied  load.  In  Figure  4.11, 
the  present  solution  is  compared  to  an  approximate  solution  given  by  Timoshenko 


[80] 


0.217 


Pa^ 

Eh* 


w, 


Wr 


=  ^+0.443(-e:) 


(4.31) 


where  Wo  is  the  transverse  deflection  at  the  center  of  a  thin  plate.  The  thickness 
and  the  radius  of  the  plate  were  0.1  in  and  2  in,  respectively.  The  Young’s  modulus 
was  3  msi  and  the  Poisson’s  ratio  was  0.3.  It  is  shown  that  the  present  solution  Eind 
Timoshenko’s  solution  agreed  fairly  well. 

An  Aluminium  Thick  Plate  Loaded  by  a  Spherical  Indenter 

The  analytical  solution  for  the  contact  force  /  with  respect  to  the  indentation 
a  is  deflned  by  the  elastic  contact  law  [10,  11]  given  by: 


/  =  ka"  (4.32) 

where  the  exponent  n  can  be  taken  as  1.5  for  a  homogeneous  transversely  isotropic 
plate  and  the  contact  stiffness  k  can  be  given  as 


k  = 


E, 

1-1/2 


(4.33) 


where  the  indenter  is  assumed  to  be  rigid.  The  i/  is  the  poisson  ration  for  the  plate 
and  the  Ez  is  the  Young’s  modulus  of  the  plate  in  the  thickness  direction.  The 
Young’s  modulus  was  3  msi  and  the  Poisson’s  ratio  was  0.3.  The  R  is  the  radius  of 
the  indenter  and  its  value  here  was  0.25  in. 
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APPLIED  DISPLACEMENT  wo  (in) 


Figure  4.11  A  load-displacement  curve  for  a  circular  plate  subjected  to  a  central 
point  load.  Comparison  between  the  cedculations  and  elzisticity  solu¬ 


tions. 
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The  plate  had  a  length  of  4  in,  a  width  of  3  in,  and  a  depth  of  1.3  in.  The  bottom 
of  the  plate  was  placed  on  the  rigid  support  (i.e.,  w  =  0  for  the  bottom  surface 
of  the  plate).  The  four  edges  were  simply-supported.  The  indenter  was  applied 
downward  at  the  center  of  the  plate.  Only  1/4  plate  was  modeled.  Figure  4.12 
gives  the  comparison  between  the  elasticity  solution  and  the  solution  obtained  by 
the  proposed  firxite  element  analysis.  It  is  shown  that  the  present  solution  and  the 
elasticity  solution  agreed  fairly  well. 
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Figure  4.12  The  load- indentation  curve  for  an  aluminium  plate  subjected  to  a 
spherical  indenter  with  a  radius  of  0.25  in.  Comparison  between  the 
predictions  and  the  analytical  solution.. 
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§4.5  Results  and  Comparison 

In  order  to  further  substantiate  the  proposed  analysis,  numerical  calculations 
were  performed  to  compare  with  the  experimentjil  data  generated  during  the  investi¬ 
gation  or  obtained  from  the  literature.  Two  types  of  material  systems  were  selected 
for  comparison:  T300/976  and  T800H/3900-2  graphite/epoxy  materials.  Both  com¬ 
posites  contain  thermoset  resin  systems,  but  T800H/3900-2  composite  is  toughened 
on  the  interfaces  with  thin  polyamide  particles  [84].  The  stiffness  and  strengths  of 
unidirectional  composites  of  both  materials  are  very  similar,  except  strain  energy' 
release  rates  Gic  and  Guc-  Due  to  its  toughening  interface,  T800H/3900-2  com¬ 
posite  has  much  higher  fracture  toughnesses  than  T300/976,  especially  for  Guc 
values  as  shown  in  Table  4-2. 

Test  data  on  T300/976  composite  were  taken  from  Firm  and  Springer  [1],  while 
experiments  on  T800H/3900-2  composite  were  conducted  during  this  investigation. 
Three  different  layups  were  selected:  [06/902]a,  [02/906],,  and  [O4/9O4],.  The  ex¬ 
periments  on  T800H/3900-2  composite  followed  the  same  procedures  as  were  used 
by  Finn  and  Springer  for  T300/976  composite.  All  the  specimens  had  the  same 
number  of  plies  (16  plies),  width  (3in),  and  length  (4in)  During  the  tests,  the  spec¬ 
imens  were  cleimped  on  the  two  parallel  edges  and  free  on  the  others.  A  transverse 
concentrated  load  induced  through  an  indenter  of  radius  0.25  in  weis  applied  by  an 
MTS  machine.  Displacement  control  was  used  throughout  all  the  tests. 

All  the  specimens  were  X-rayed  before  and  after  each  test  to  examine  internal 
damage  resulting  from  the  applied  load.  The  applied  energy  (load)  as  a  fimction 
of  measured  delamination  size  from  X-radiographs  was  also  recorded.  Figure  4.13 
shows  typical  damage  patterns  from  X-radiographs  for  three  different  layups  of 
T300/976  composite  at  a  given  input  energy.  The  straight  lines  indicated  matrix 
cracks  and  the  area  bounded  by  white  contour  curve  was  delsunination.  It  is  worth 
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noting  that  the  shape  of  the  delamination  strongly  depended  upon  the  ply  orienta¬ 
tion  and  was  considerably  different  for  each  of  the  three  layups. 


Table  4-2 

Material  properties  for  Fiberite  T800H/3900-2. 


Material  Property 

Value 

Units 

Longitudinal  Yoxmg’s  modulus,  Ex 

23.2 

msi 

Trems verse  Young’s  modulus,  Ey 

1.33 

msi 

Out-of-plane  Yoimg’s  modulus.  Ex 

1.33 

msi 

Poisson’s  ratio,  Vxy 

0.28 

Poisson’s  ratio,  Vxz 

0.28 

Poisson’s  ratio,  Vyx 

0.28 

In-plane  shear  modulus,  Gxy 

0.90 

msi 

Out-of-plane  shear  modulus,  G** 

0.90 

msi 

Out-of-plane  shear  modulus,  Gy* 

0.90 

msi 

Longitudinal  tensile  strength,  Xt 

413 

ksi 

Longitudinal  compressive  strength,  Xc 

191 

ksi 

Transverse  tensile  strength,  Yj 

6.41 

ksi 

Transverse  compressive  strength,  Yc 

24.4 

ksi 

LongitudinaJ  shear  strength  (cross  ply),  Sep 

52.89 

ksi 

Critical  strain  energy  release  rate  for  mode  I,  Gj^ 

1.50 

Ibf/in 

Critical  strain  energy  release  rate  for  mode  II,  G/Zp 

18.0 

Ibf/in 

Ply  thickness  hg 

0.006535 

in 
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[06/902]s 


[02/906]s 


[04/904]s 


Figure  4. 


T300/976 


X- Radiographs  of  graphite/epoxy  composites  resulting  from  a  quasi- 
static  transverse  load,  a)  [Og/QOa)#  layup;  b)  [02/906]*;  and  c)  [O4/9O4]*. 
Data  taken  from  [1]. 
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For  [06/902]j  layup,  a  long  slender  peanut-like  shape  of  a  delzunination  on  the 
interface  near  the  bottom  surface  was  observed.  A  matrix  crack  (bending  crack) 
in  the  bottom  0  degree  ply  group  was  also  clearly  seen  from  the  X-radiograph 
given  in  Figure  4.13.  The  left  and  right  half  of  the  delaminations  always  grew 
in  a  relatively  equal  size.  However,  for  [02/906]*  layup,  the  damage  pattern  was 
considerably  different  from  [06/902]*.  Distinct  matrix  cracks  (shear  cracks)  in  the 
centred  90  degree  ply  group  near  the  loading  area  could  be  seen  cleeiry  in  addition 
to  the  bending  crack,  and  the  delamination  shape  was  much  wider  and  close  to  a 
semi-circular  shape  on  each  side  which  was  quite  unevenly  distributed.  Meeinwhile, 
the  [O4/9O4]*  layup  seems  to  produce  a  damage  pattern  somewhere  in  between  the 
[O6/9O2]*  and  the  [O6/9O2]*  layups. 

For  T800H/3900-2  composite,  the  damage  patterns  were  very  similar  to  those 
of  T300/976  composite  for  the  same  ply  orientation,  except  that  the  deleimination 
sizes  were  substantially  smaller  than  those  of  T300/976.  There  was  no  delamination 
damage  below  a  certain  energy  level  in  T800H/3900-2  composite,  at  which  the 
T300/976  composite  weis  already  severely  damaged.  At  a  given  amount  of  the 
applied  load,  when  panels  made  of  both  materials  had  damage,  the  damage  size  in 
[O2/9O6]*  of  T800H/3900-2  was  about  only  one- fifth  of  the  observed  from  T300/976. 
A  comparison  of  damage  size  between  T300/976  and  T800H/3900-2  for  various 
amounts  of  input  energy  is  shown  in  Figure  4.14. 

TSOO/976  [O6/9O2]*  Panels 

In  order  to  evaluate  the  damage  propagation  in  [O6/9O2]*  composites  due  to 
transverse  concentrated  loading,  both  types  1  and  2  damage  models  were  utilized 
to  perform  the  numerical  simulations.  The  type  1  model  considered  only  a  de- 
Izunination  without  a  bending  crack,  however,  type  2  model  considered  both  the 
delamination  and  the  bending  crack.  The  typical  meshes  used  in  the  calculations 


DELAMINATION  SIZE 


Chapter  4:  Damage  Induced  by  a  Spherical  Indenter 


96 


Figure  4.14  Comparison  between  the  measured  and  predicted  delamination  sizes 
in  [Oe/OOa],  composites  made  of  T300/976  and  T800H/3900-2  with 
respect  to  quasi-impact  energy. 
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were  given  previously  in  Figure  4.13-a. 

The  mesh,  deformation,  and  stress  data  generated  from  the  finite  element  anal¬ 
ysis  were  stored  in  standard  post-processing  files  and  transformed  to  the  commercial 
code  I  -  DBAS  [69]  (licensed  by  Structural  Dynamics  Research  Corporation).  Fig¬ 
ure  4.15  shows  a  close-up  view  of  a  deformed  cross-section  {X  —  Z,  Y  =  0)  of  a 
[06/902]«  composite  at  a  load  of  172  Ibf.  The  local  indentation  due  to  the  rigid 
indenter  can  be  clearly  observed  in  the  figure.  A  relatively  large  sliding  between  the 
upper  and  lower  surfaces  of  the  delamination  were  also  clearly  shown  in  the  figure. 

A  close-up  side  view  of  a  deformed  cross-section  (Y  —  Z,  X  =  0)  of  a.  [06/902]j 
is  shown  in  Figure  4.16.  The  bending  crack  was  open  due  to  the  applied  load,  which 
closely  resembles  the  cylindrical  bending  of  a  [90m/0n]»  beam  containing  a  surface 
crack.  The  local  indentation  due  to  the  rigid  indenter  and  relatively  large  sliding 
between  the  upper  and  lower  surfaces  of  the  delamination  were  also  shown  in  the 
figiire. 


[06/902].  composite 
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Figure  4.17  shows  the  comparison  of  the  calculated  total  strain  energy  release 
rate  Gt  based  on  types  1  and  2  models  for  a  laminate  containing  a  smEdl  circular 
delamination  with  and  without  the  surface  matrix  crack  (bending  crack).  For  the 
laminate  containing  the  matrix  crack,  the  value  of  along  the  delamination  front 
near  the  location  of  the  the  matrix  crack  increased  significantly  and  reached  a 
peak  at  the  intersection  between  the  matrix  crack  and  delamination  (at  =  0“). 
The  value  of  Gt,  however,  decreased  rapidly  as  the  delamination  front,  where  the 
value  of  Gt  was  calculated,  moved  away  from  the  location  of  the  matrix  crack  {<f) 
approaches  to  90°). 

A  completely  different  distribution  of  the  total  strain  energy  release  rate  was 
obtained  for  the  laminate  containing  no  bending  crack.  Overall,  the  values  of  Gr 
obtained  from  the  type  1  model  were  much  smaller  them  those  calculated  by  the 
type  2  model.  Apparently,  the  laminate  with  the  matrix  crack  could  initiate  the 
delamination  growth  at  a  much  earlier  loading  stage  than  the  one  without.  Once 
the  delamination  propagated,  it  would  grow  along  the  direction  of  matrix  cracking. 

Figure  4.18  shows  the  sequence  of  the  deleimination  growth  in  the  laminate 
without  the  matrix  crack  predicted  by  the  type  1  model  based  on  the  crack  growth 
criterion  (Eq.  (4.8)).  The  number  shown  in  the  bracket  at  the  upper  right  corner 
of  each  sub-figure  indicates  the  sequence.  The  sub-figures  were  generated  by  eval¬ 
uating  the  strain  energy  release  rate  ratio  Ej  for  various  sizes  of  delaminations. 
First,  the  values  of  Ed  were  calculated  for  a  small  circular  implanted  delamination. 
If  delamination  growth  was  predicted  {Ed  >1),  then  the  delamination  size  was 
extended  slightly  along  its  major  (X-direction,  <l>  =  0°)  or  minor  (F-direction,  <j)  = 
90°)  axis  according  to  the  predicted  direction  of  delamination  growth.  In  this  anal¬ 
ysis,  the  delamination  was  only  adlowed  to  grow  into  either  a  circular  or  an  elliptical 
shape.  Again,  numerical  calculations  were  re-performed  to  evaluate  the  strain  en- 
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Figure  4.17  Comparison  of  the  calculated  total  strain  energy  release  rate  along  a 
delamination  front  in  a  [06/902]*  composite  with  or  without  a  surface 
matrix  crack. 
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ergy  release  rates  for  the  laminate  with  the  new  delamination.  This  procedure  wzis 
repeated  until  the  calculated  strain  energy  release  rate  ratio  was  smedler  than  unity 
everywhere  along  the  delamination  front.  The  size  of  the  final  delamination  wzis 
then  considered  to  be  the  delamination  shape  corresponding  to  the  given  loading 
condition. 

At  0.05  inches  of  indenter  displacement,  Figure  4.18-1  shows  that  delamination 
growth  initiated  at  ^  =  0®.  As  the  delamination  expanded  along  its  major  axis  from 
Figure  4.18-1  to  Figure  4.18-2,  the  strain  energy  release  rate  increased,  indicating 
2in  unstable  growth.  It  is  noted  that  as  the  delamination  continued  to  expand  along 
its  major  axis,  the  peak  of  the  strain  energy  release  rate  ratio  gradually  shifted 
more  and  more  toward  =  90°,  indicating  a  much  more  imiform  expansion  of 
delamination.  Figure  4.18-3  indicates  delamination  expansion  along  the  minor  axis 
due  to  high  strain  energy  release  rates  near  =  90°.  The  distribution  of  the  strain 
energy  release  rate  ratio  along  the  delamination  front  for  all  the  delaminations  was 
relatively  uniform  compared  to  that  of  the  Isiminate  with  the  matrix  crack.  As  a 
result,  the  delamination  tended  to  grow  into  an  elliptical  shape  with  major  and 
minor  axis  ratio  about  2. 

However,  the  deleimination  shape  predicted  by  the  type  1  model  was  signifi¬ 
cantly  different  from  the  shape  observed  in  the  experiment.  First,  the  well-known 
peanut-like  delamination  could  not  be  predicted.  Second,  the  predicted  shape  was 
fairly  large  instead  of  slender  and  was  very  different  from  the  observed  delamination 
shape  with  major  and  minor  axis  ratio  about  5.  And  third,  the  applied  load  (or 
displacement)  which  initiated  the  delamination  growth  was  much  higher  than  the 
experimental  one.  It  was  noted  that  something  was  missing  in  the  type  1  model  for 
simtilating  the  delamination  growth  due  to  a  spherical  indenter  loading.  It  seemed 
to  show  that  the  type  1  model  could  not  predict  accurately  the  damage  growth  in 
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Figure  4.18  The  predicted  delamination  growth  sequence  in  a  [06/902]j  composite 


without  a  surface  matrix  crack. 
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this  laminate  due  to  a  spherical  indenter. 

However,  a  completely  different  growth  pattern  was  predicted  by  the  type  2 
model.  First,  at  the  same  load,  delamination  growth  was  predicted  to  occur  much 
earlier  for  the  laminate  containing  the  same  small  delamination  given  in  Figure  4.19- 
1.  Figure  4.19  shows  the  sequence  of  the  delamination  growth  in  the  laminate  with 
the  matrix  crack  predicted  by  the  type  2  model. 

At  a  fixed  indenter  displacement.  Figure  4.19-1  and  Figure  4.19-2  show  that 
the  initial  delamination  grew  from  a  small  circular  shape  into  a  slender  elliptical 
shape  along  the  direction  at  <^  =  0°  with  its  major  axis  parallel  to  the  fiber  direction 
of  the  bottom  ply  group.  The  delamination  growth  was  unstable  because  the  strain 
energy  release  rates  actually  increased  as  the  delamination  expanded.  Figure  4.19- 

3  indicates  that  delamination  continued  to  expand  along  its  major  axis,  until  the 
strain  energy  release  rate  ratio  stzirted  to  decreaise  along  the  delamination  front  near 
4>  =  0°. 

After  a  substantial  growth  of  the  delamination  along  its  major  axis.  Figure  4.19- 

4  shows  that  the  strain  energy  release  rates  started  to  increase  near  <f>  =  75°  away 
from  the  matrix  crack,  initiating  the  expansion  of  the  delamination  along  the  minor 
axis.  However,  a  slight  expansion  of  the  delamination  along  its  minor  axis  as  shown 
in  Figure  4.19-5  subsequently  caused  a  significant  increase  of  the  strain  energy 
rele^lse  rates  along  the  major  axis.  Hence,  the  deleunination  started  to  grow  along  the 
major  axis  again.  Apparently,  the  growth  of  the  delamination  was  predominantly 
controlled  by  the  delamination  front  near  the  neighborhood  where  the  matrix  crack 
intersected  with  the  delamination.  It  is  worth  noting  that  at  <j>  =90°,  strain  energy 
release  rates  remained  at  a  minimum  regardless  of  the  : '  '.pe  of  the  delamination. 
This  observation  implies  that  no  delamination  would  grow  in  this  area,  leading  the 
delamination  to  a  peanut  shape  which  was  frequently  observed  from  the  experiments 
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Figure  4.19  The  predicted  delamination  growth  sequence  in  a  [06/902)5  composite 
containing  a  surface  matrix  crack. 
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for  this  type  of  ply  orientataion,  such  as  the  one  given  in  Figure  4.13. 

Figure  4.20  shows  the  comparisons  between  the  estimated  delamination  sizes  in 
the  X  and  y  directions  emd  the  sizes  measured  from  the  experiments.  The  correlation 
between  the  experiments  and  the  predictions  was  very  good.  Figure  4.21  shows 
the  calculated  strain  energy  release  rates  of  Modes  I,  II,  and  III  for  the  laminate 
contziining  a  small  and  a  large  delamination  with  a  surface  matrix  crack.  Cleaxly, 
for  the  laminate  containing  the  initial  surface  matrix  crack,  Gj  (Mode  I  fractme) 
dominated  the  total  strain  energy  release  rate  along  the  delamination  front  near  the 
neighborhood  where  the  surface  matrix  crack  intersected  with  the  delamination  {6 
=  0°).  This  observation  is  also  valid  for  a  moderately  large  delamination,  as  shown 
in  the  same  figure. 

However,  when  the  surface  crack  was  ignored  in  the  analysis.  Modes  II  and 
III  fractures  dominated  the  total  strain  energy  release  rate  along  the  delamination 
front  for  the  laminate,  as  shown  in  Figure  4.22.  Although  the  contribution  of  each 
mode  to  the  growth  of  the  delamination  strongly  depended  upon  the  current  shape 
of  delamination,  the  presence  of  the  matrix  crack  cleaxly  played  a  very  important 
role  in  the  delamination  growth.  Therefore,  it  is  very  importzint  that  the  initial 
matrix  crack  be  considered  in  the  analysis  for  understamding  the  damage  mechanics 
and  mechanism  of  laminated  composites  due  to  transverse  loading. 

T300/976  [O2/9O6],  Panels 

Type  3  and  4  models  were  used  to  study  the  delamination  growth  behavior 
in  [02/906]j  composites.  The  type  3  model  simulates  a  delamination  induced  by  a 
pair  of  shear  cracks,  while  the  type  4  model  considers  both  a  bending  crack  and 
shear  cracks  in  conjimction  with  a  delamination.  Therefore,  numerical  calculations 
were  generated  for  the  laminate  conteuning  a  deleimination  with  veuious  sizes  and  a 
shear  crack  in  the  middle  90®  ply  group  with  and  without  a  surface  bending  crack 
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Figure  4.21  Compzurison  of  the  calculated  strEun  energy  releeise  rates  of  Modes  I, 


II,  and  III  along  a  small  and  a  laurge  delamination  front  in  a  [06/902]* 
composite. 
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in  the  bottom  0®  ply  group.  As  suggested  by  the  experimental  observation,  the 
lengths  of  the  matrix  cracks  were  assumed  to  be  slightly  larger  than  the  delamination 
length  and  width  in  the  i  and  y  directions  throughout  the  calculations.  In  the 
following  figures,  similar  to  the  [Oe/OOa]*  case,  numericsil  results  will  be  presented 
for  delaminations  from  small  to  large  sizes. 

A  close-up  view  of  a  deformed  cross-section  {X  —  Z,Y  =  0)  of  a  [02/906]a  panel 
at  a  load  of  327  Ibf  is  shown  in  Figure  4.23.  The  local  indentation  due  to  the  rigid 
indenter  eind  relatively  large  sliding  between  the  two  surfaces  of  the  delamination 
and  internal  crack  are  clearly  shown  in  Figure  4.23. 

Figure  4.24  shows  the  sequence  of  the  delamination  growth  in  the  Izuninate 
contziining  a  pair  of  shear  cracks  predicted  by  the  type  3  model.  At  the  fixed  indenter 
displacement,  Figvire  4.24-1  and  Figure  4.24-2  show  that  the  initial  delzimination 
grew  from  a  small  elliptical  shape  into  a  wider  elliptical  shape  eilong  the  direction 
at  </»  =  90°  with  its  minor  axis  perpendicular  to  the  fiber  direction  of  the  bottom 
ply  group.  The  delamination  growth  was  unstable  because  the  strain  energy  release 
rates  increased  as  the  delamination  expanded.  Figure  4.24-3  indicates  that  the 
delamination  begain  to  expeind  along  its  major  axis,  as  the  strain  energy  release  rate 
ratio  started  to  decrease  aJong  the  delamination  front  ne£ir  (f)  =  90°  and  had  the 
highest  value  at  the  delamination  front  near  ^  =  0°. 

The  expansion  of  the  delamination  along  its  major  or  minor  eocis  ceased  to 
increase  when  the  growth  index  Ed  along  the  delamination  front  was  below  unity,  as 
in  the  case  shown  in  Figure  4.24-6.  The  delamination  shape  shown  in  Figure  4.24- 
6  is  approximately  a  circular  one,  which  agreed  with  the  delamination  shape  of 
[02/906]*  panels  observed  from  the  experiments  (see  Figure  4.13-b).  Figure  4.25 
shows  the  comparisons  between  the  estimated  delamination  sizes  in  the  x  and  y 
directions  and  the  sizes  measured  from  the  experiments.  The  correlation  between 
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Figure  4.24  The  predicted  delaxnination  growth  sequence  in  a  [02/906],  composite 


containing  a  surface  matrix  crack  and  an  internal  matrix  crack. 
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the  experiments  and  the  predictions  was  good. 

It  was  interesting  to  evaluate  the  contributing  fractmre  modes  for  the  delami¬ 
nation  growth  in  [02/906]*  panels.  Figure  4.26  shows  the  calctilated  strain  energy 
release  rates  of  Modes  I,  II,  and  III  for  the  laminate  containing  a  small  delamination 
with  an  internal  shear  crack  and  a  surface  bending  crack.  For  the  laminate  with 
the  initial  bending  crack,  G /  (Mode  I  fractme)  contributed  significantly  to  the  total 
strain  energy  release  rate  along  the  delamination  front  near  the  neighborhood  where 
the  matrix  cracks  intersected  with  the  delamination  {<f>  =  0“  and  90“),  indicating 
that  delamination  initiation  was  in  mixed  modes. 

Figme  4.27  shows  the  calculated  strain  energy  relezise  rates  of  Modes  I,  II, 
and  III  for  the  laminate  containing  a  larger  delamination  with  an  internal  crack 
and  with  or  without  a  surface  matrix  crack.  Clearly,  Modes  II  and  III  fractures 
dominated  the  total  strain  energy  release  rate  for  the  laminate  with  and  without 
the  surface  matrix  crack.  And  it  is  noted  that  the  total  strain  energy  release  rate 
was  only  slightly  affected  by  the  introduction  of  the  surface  crack,  indicating  that 
delamination  growth  was  donoinated  by  the  interned  shear  crack. 
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Figure  4.26  Comparison  of  the  calculated  streiin  energy  release  rates  of  Modes  I, 
II,  and  III  along  a  smadl  delamination  front  in  a  [02/906]*  composite 
containing  am  internal  crack  amd  a  surface  matrix  crack. 
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Figure  4.27  Comparison  of  the  calculated  strain  energy  release  rates  of  Modes  I, 
II,  and  III  along  a  large  delamination  front  in  a  [Oa/QOe]*  composite 
containing  an  internal  crack  and  with  or  without  a  surface  matrix 
creick. 
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TSOO/976  [O4/9O4],  Panels 

It  hzis  been  shown  that  delamination  growth  in  [06/902)5  peniels  was  meiinly 
dominated  by  Mode  I  fracture,  while  in  [Oa/OOe]*  it  was  mainly  dominated  by  shear 
(Mode  II  and  III)  fracture.  In  the  following,  delamination  growth  of  a  [04/904)5 
panel  was  studied.  Again,  type  3  and  4  models  were  considered. 

Figure  4.28  shows  the  sequence  of  the  delamination  growth  in  the  laminate  with 
the  matrix  cracks  predicted  by  the  type  3  model.  At  a  fixed  indenter  displacement. 
Figure  4.28-1  and  Figure  4.28-2  show  that  the  initial  delamination  grew  from  a  small 
elliptical  shape  into  a  larger  elliptical  shape  along  the  direction  at  (^  =  90°  with  its 
minor  axis  perpendicular  to  the  fiber  direction  of  the  bottom  ply  group. Figure  4.28-3 
indicates  the  delamination  would  begin  to  expand  along  its  major  axis,  eis  the  strain 
energy  release  rate  ratio  started  to  decrease  along  the  delamination  front  around  4> 
=  90°,  and  had  the  highest  values  at  the  delamination  front  around  </►  =  0°. 

Similarly,  the  deleimination  ceased  to  grow  along  its  major  or  minor  axis  as  the 
growth  index  Ed  along  the  delEimination  front  fell  to  below  unity,  such  as  the  case 
shown  in  Figure  4.28-4.  The  delamination  shape  shown  in  Figure  4.28-4  wets  ein 
elliptical  shape,  which  agreed  with  the  shape  observed  in  the  experiments  (see  Fig¬ 
ure  4.13-c).  Figure  4.29  shows  the  comparisons  between  the  estimated  delamination 
sizes  in  the  x  and  y  directions  and  the  sizes  measured  from  the  experiments.  The 
correlation  matches  reasonably  well  between  the  experiments  and  the  predictions. 

It  was  again  interesting  to  investigate  the  contributing  modes  to  the  delami¬ 
nation  growth.  Figure  4.30  shows  the  comparisons  of  the  calculated  crack  growth 
index  Ed  for  the  laminate  containing  a  small  elliptical  delzunination  with  em  internal 
matrix  crack  and  with  or  without  the  surface  matrix  crack.  The  upper  figure  was 
the  case  for  a  small  delamination.  For  the  laminate  containing  the  surface  matrix 
crack,  the  value  of  Ed  along  the  delamination  front  near  the  location  of  both  matrix 
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Figure  4.28  The  predicted  delamination  growth  sequence  in  a  [04/904]*  composite 
contzdning  a  surface  matrix  crack  emd  an  interned  matrix  crack. 
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Figure  4.29  The  comparison  of  delamination  sizes  between  the  predicted  and  the 
measured  in  a  [O4/9O4]*  composite  under  a  spherical  indenter. 
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cracks  increased  significantly  and  reached  a  peak  at  the  intersections  between  the 
matrix  cracks  and  delamination  (at  4>  —  90°).  The  value  of  Ei,  however, 

decreased  rapidly  as  the  delamination  front  where  the  value  of  Ed  was  calculated 
moved  farther  away  from  the  location  of  the  matrix  cracks  (when  <l>  is  not  around 
0°  and  90°).  Therefore,  for  a  small  delamination,  the  coupling  between  the  matrix 
cracks  and  delamination  was  significant. 

A  different  distribution  of  the  total  strain  energy  release  rate  was  obtained 
for  the  laminate  containing  no  pre-introduced  matrix  crack.  Overall,  the  values 
of  Ed  were  much  smaller  than  those  calculated  from  the  laminate  containing  the 
matrix  crack.  Apparently,  the  laminate  with  the  surface  matrix  crack  could  initiate 
the  delamination  at  a  much  earlier  loading  stage  than  the  one  without.  Once  the 
delamination  propagated,  it  would  grow  along  with  the  direction  of  internal  matrix 
cracking. 

However,  for  moderate  and  large  delaminations  such  as  the  case  shown  in  the 
lower  figure  of  Figure  4.30,  the  surface  crack  had  limited  contribution  to  the  delami¬ 
nation  growth.  The  coupling  between  the  surface  matrix  crack  and  the  delamination 
only  appeeired  locally  at  <^  =  0°.  Apparently,  the  effect  of  the  surface  matrix  crack 
on  the  delamination  growth  decreased  significantly  eis  the  delamination  grew. 

It  was  also  interesting  to  examine  the  contributing  fracture  modes  for  delam¬ 
ination  growth  when  the  delamination  is  small  or  moderately  sized.  Figure  4.31 
shows  the  ceilculated  strain  energy  release  rates  of  Modes  I,  II,  and  III  for  the  lam¬ 
inate  containing  a  small  delamination  with  an  internal  crack  and  a  surface  matrix 
crack.  Clearly,  for  the  laminate  with  a  small  delamination,  Gj  (Mode  I  fracture) 
contributed  significantly  to  the  totzd  strain  energy  release  rate  along  the  delamina¬ 
tion  front  near  the  neighborhood  where  the  surface  matrix  crack  intersected  with 
the  delamination  (<^  =  0°  and  90°).  Genereilly,  the  growth  for  a  small  delamination 
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Figure  4.30  The  effect  of  bending  cr£w:k  on  the  crack  growth  index  for  a  laminate 
containing  a  delamination  eind  a  bending  and  shear  crack.  Comparison 
of  Ed  vzdues  between  a  small  and  a  large  delamination  in  a  [O4/9O4], 
composite. 
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Figure  4.31  The  calculated  strain  energy  release  rates  between  a  small  and  a  large 
delamination  front  in  a  [04/904]s  composite  containing  an  internal 
crack  and  a  surface  matrix  crack. 
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was  in  mixed  modes.  However,  the  lower  part  of  Figure  4.31  showed  that  except 
in  the  neighborhood  of  ^  =  0®,  the  strain  energy  release  rates  were  only  slightly 
affected  by  the  surface  crack.  The  shear  mode  fractures  (Modes  II  and  III)  took 
the  consideration  of  the  largest  greatest  portion  of  the  total  energy  release  rate,  in¬ 
dicating  that  the  delamination  growth  was  dominated  by  interned  matrix  crewdcing 
and  the  surface  crack  effect  was  minimmn. 

T800H/S900-2  [Oj/OOe],  Panels 

Numerical  simulations  were  also  performed  on  T800H/3500-2  composites  to 
estimate  delamination  sizes  corresponding  to  the  applied  load  (energy).  The  results 
of  the  test  data  on  [Oa/OOe]*  are  shown  in  Figure  4.14.  Only  the  type  4  model  was 
used  to  correlate  the  predictions  and  the  experiments.  The  estimated  delamination 
sizes  in  the  x  and  y  directions  are  presented  in  Figure  4.14.  The  comparison  between 
the  estimated  sizes  and  the  measured  sizes  agree  quite  well. 

For  the  purpose  of  comparison,  the  estimated  results  and  experimental  data 
for  T300/976  panels  are  also  shown  in  the  same  figure.  It  is  shown  clearly  that  the 
panels  made  of  T800/976  have  much  smaller  damage  sizes  than  the  panels  made  of 
T300/976.  It  was  also  observed  that  the  damage  shape  was  much  more  symmetric 
than  the  that  in  T300/976  panels,  indicting  the  damage  growth  inside  the  toughened 
composites  tend  to  have  more  stable  growth.  The  reduced  damage  size  is  one  of 
the  main  reasons  that  contribute  to  the  reported  high  compressive  strength  after 
impact  (CAI)  for  this  new  material  system. 
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A  finite  element  analysis  was  developed  for  analyzing  cross-ply  laminated  com¬ 
posite  plates  subjected  to  transverse  concentrated  loading  resulting  from  a  spherical 
indenter.  Based  on  the  analysis,  the  following  remarks  ctin  be  made  for  the  lami¬ 
nates  studied: 

1.  The  initial  matrix  cracking  significantly  affects  the  growth  of  delamination 
resulting  from  transverse  loading. 

2.  Matrix  cracks  enhance  delamination  growth. 

3.  Mode  I  fracture  is  critical  for  surface  crack-induced  delamination  growth. 

4.  Modes  II  and  III  shear  fracture  dominate  internal  crack  induced  delamination 
growth. 

5.  Matrix  cracking  must  be  considered  for  modeling  delamination  propagation 
due  to  transverse  loading. 

6.  Delamination  growth  depends  on  ply  orientation. 

Additionally,  it  is  noted  that  the  proposed  finite  element  analysis  cein  be 
extended  to  study  other  ply  orientations  and  other  loading  conditions,  such  as 
delamination-buckling. 
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A  2-D  Finite  Element  Mesh  Remodeling  Procedure 


Before  a  structure  is  analyzed  or  tested,  there  is  no  information  about  where 
the  matrix  cracks  are  to  occur  or  what  interfaces  are  to  be  delaminated.  Here  a 
simple  but  effective  scheme  is  developed,  which  can  be  illustrated  by  two  simple  but 
typical  meshes,  shown  in  Figures  A.l  and  A.2.  The  first  mesh,  shown  in  Figure  A.l, 
is  mainly  used  for  [0m/90n]«  laminates,  for  which  embedded  crack-induced  delam¬ 
inations  are  the  typical  damage  mode.  The  second  mesh,  shown  in  Figure  A. 2, 
is  mainly  used  for  [90m/0n]s  laminates,  for  which  surface  bending  crack-induced 
delzunination  is  the  typical  damage  pattern. 

This  technique  uses  the  originzd  two-dimensional  finite  element  grid  system.  In 
the  original  mesh,  double  lines  instead  of  one  regular  line  are  iised  on  eaeh  interface 
between  elements  in  the  thickness  direction  or  in  the  span-length  direction.  In 
other  words,  each  node  is  assigned  systematicedly  four  nodal  numbers,  instead  of 
one,  in  the  global  nodal  system.  If  no  damage  occurs,  adl  four  numbers  at  each 
node  are  assigned  the  same  equation  numbers.  For  example,  the  four  nodes  i,j,  k,  I 
in  Figure  A.l  have  the  same  equation  numbers  if  there  is  no  damage  introduced  at 
that  point.  Once  a  crack  is  generated,  nodes  on  different  sides  of  the  double  lines 
are  split  and  therefore  have  different  equation  numbers.  For  example,  eis  shown  in 
Figure  A-1,  nodes  {E,  G)  and  nodes  {F,  H)  Eire  separated  due  to  the  existence  of  the 
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vertical  matrix  crack,  while  two  nodes  E  and  G  (same  for  F  smd  H)  are  assigned 
the  same  equation  numbers.  For  some  corner  points  such  as  (A,  B,  C,  B),  where 
the  delamination  and  the  matrix  crack  intersect,  (A,B,D)  are  zissigned  the  same 
equation  number.  C  remains  independent. 

As  the  matrix  cracks  and  delaminations  may  intersect,  one  needs  to  apply  the 
contact  condition  to  the  intersection  points  (such  as  C  and  J)  to  prevent  the  over¬ 
lapping  of  the  matrix  crack  and  the  delamination  surfaces.  Therefore,  five  degrees 
of  freedom  need  to  be  assigned  to  each  node.  The  first  three  are  for  displacements  in 
Xi,  X2,  and  directions.  The  fourth  is  used  for  a  Lagrzmge  multiplier  for  possible 
delaminations  and  the  fifth  is  used  for  a  Lagrange  multiplier  for  possible  matrix 
cracks. 

By  this  technique,  matrix  cracking  can  be  easily  simulated  wherever  it  occurs, 
as  the  interface  only  needs  to  be  split  between  two  elements  and  the  crack  is  intro¬ 
duced  without  much  difficulty.  This  is  also  true  for  delamination  modeling. 
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Figure  A. 2  A  typical  mesh  used  in  the  analysis  for  modeling  bending  matrix  crack 
induced  delamination  in  two  dimensions. 
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Contact  Residual  and  Stiffness 
of  a  2-D  Contact  Element 


In  order  to  assure  quadratic  convergence  of  Newton’s  method,  it  is  essential 
to  derive  a  consistent  tangent  matrix  in  lineeirizing  the  nonlinear  finite  element 
equation.  In  the  following,  a  contact  element  consisting  of  one  slave  node  and  two 
master  nodes  is  derived.  The  derivation  is  in  a  vector  and  tensor  form.  In  order  to 
make  the  implementation  easier,  matrix  expressions  are  derived  for  contact  residuzJ 
and  consistent  contact  tangent  stiffness  matrix. 

By  referring  to  Figure  B.l,  the  contact  gap  associated  with  a  typical  slave  node 
s  is  given  by  a  simple  form 


g  =  n-ix,-Xi) 


(B.l) 


where  n  is  the  normal  of  contact  point  c  on  master  segment  1  —  2.  It  is  worthwhile 
to  note  that  the  normal  is  unique  for  master  segment  1  —  2.  To  obtain  an  explicit 
expression  for  the  residual  and  contact  stiffness  in  Eq.  (4.20)  and  Eq.  (3.36),  the 
variation  6g  must  be  performed.  In  the  fully  nonlinear  ceise,  the  cheinge  in  the  normal 
vector  n  needs  to  be  considered.  As  Tl  =  62  x  ei,  the  change  in  the  vector  Cj  also 
needs  to  be  considered.  The  tzmgent  Ci  is  given  by  61  =  (*2  “  3;i)/||®2  ~  ®i||- 
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Figure  B.l  Geometry  of  a  slave  node  s  with  a  master  segment  1-2. 
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By  using  the  chain  rule,  the  variation  of  61  is  given  as 
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- [iT [(^®2  -  ^*1)11*2  -  «i||  -  (®2  -  a;i)^lla;2  -  ®i||]  {B.2) 

ll®2  -  «l|| 

Recalling  ||aJ2  —  ®i||  =  [(*2  —  ®i)  •  (®2  ~  the  following  is  obtained: 


6||a;2-a;il|  = 


{X2  -  JCi)  •  {8x2  -  ^aii) 
!l®2-a;ill 


(B.3) 


Therefore 

■||  [(^g2  -Sx^)-  {X2  -  •  {SX2  -  <5ai)]  {BA) 

®2  ~  aiii  ajj  —  ai 


=  TZ - -  ei  ei  •  {6x2  -  Sx^)] 

X2  -  ai 


(5.5) 


— zn\(^  -  ci  ®  ei){6x2  -  sxi) 

X2  —  ai  I 


(5.6) 


=  TZ - ZH^'^  ®  n){6x2  -  6x1) 

X2  -  ai 


(5.7) 


where  n  ®  Tl  +  61  ®  61  =  I  is  used.  The  symbol  ®  denotes  tensor  product.  Also, 
the  variation  of  Tl  gives  the  expression: 


6n  =  Se2  X  Cl  +  62  X  Sci  =  62  x  Sci 


(5.8) 


=  62  X  p— ^(n  ®  n)(Sx2  -  Sxi) 
Recedling  C2  x  n  =  —  Ci,  the  following  is  obtained: 


(5.9) 


— :r7r(®i  ®  n)(Sx2  -  Sxi) 

X2  —  Xjfl 


(B.IO) 
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The  variation  of  the  gap  g  therefore  can  be  derived  as  follows: 

6g  =  {8x,- Sxi)-n  +  {Xs-X\)  ■  Sn  (-B-H) 

or 

6g  =  [(6®,  -  Sxi)  •  n  -  1,^"  ~  ®  n)(6®2  -  ^*1)]  {B.12) 

||®2  ~  *lll 

or 

6g  =  [(5®,  -  6®i)  •  n  -  •  CiTi  •  (6x2  -  ^®i)]  (B.13) 

||®2  -  ®i|l 

It  is  observed  that  (®,-®i)/||®2  —  ®i||-ei  appearing  in  Eq.  (B  13)  is  exactly  the 
line  coordinate  as  can  be  seen  from  Figure  B.l.  The  contact  node  is  projected 
onto  line  1-2  and  interpolated  in  terms  of  ®i  and  ®2  as  ®c  =  (1  — 

Therefore  is  rewritten  as 

6g  =  {SXs  -  {I ^^2) -n  (B.14) 

or 

6g  =  6{Xj  -  Xc)  •  n  (B.15) 

The  total  differential  of  6g  can  be  written  as 

d(6g)  =  d[6{x,  -  ®c))  •  n  -t-  8{x,  -  ®c)  •  dn  (B.16) 

where  the  second  term  on  the  RHS  side  of  the  Eq.  (B.16)  is  simple  to  find.  The 
first  term  of  the  RHS  side  of  the  Eq.  (B.16)  is 

d[^(®,  -  ®c)]  •  n  =  -d(®c)  •  n  =  -(8Xc,(d^)  ■  n  (5.17) 

By  the  definition  of  ®c  in  terms  of  ®i  and  ®2  (that  is,  ®c  =  (1  —  +  ^*2)1  the 

following  is  obtained: 


8Xc,(  =  8x2-8Xi 


(B.\S) 
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The  differential  of  ^  =  (®,  —  Xi)/||X2  —  Xi||  •  is  derived  by  chain  rule  as 


d^  = 


1 


i®2  -®ll 
+  dCi  •  (du,  -  dUc)] 


■{gn  ■  (dU2  -  dUi) 


(B.19) 


With  all  the  information  provided  above,  the  total  differential  of  6g  can  be 
derived  as 


diSg)  = 


-  tic)  •  (ei  (8)  n){dU2  -  dui) 


\X2  -  ajii 
-  6(ii,  -  tic)  •  (Ci  (gi  n)(dti2  -  dui ) 
9 


+ 


\X2-Xi 


:d{U2  -  til)  •  (tt  ®  n)(^tl,  -  8Uc) 


(B.20) 


The  consistent  contact  stiffness  then  follows  from  the  identity 


dlSg)  =  du^  ^  ^  SU^ 

^  du^du^ 


(5.21) 


From  the  virtual  work  linearization,  the  contact  stiffness  can  be  written  as 


k  -  A— 

^  ~  du^ 


) 


And  the  contact  residual  force  of  the  element  is 


fc  = 


A 


du^ 


iB.22) 


(B.23) 


In  the  following,  explicit  matrix  expressions  r»rr  derived  for  a  single  slave  node 
in  contact  with  the  master  segment.  Based  on  these  matrices,  a  standard  assembly 
procedure  can  be  used  to  add  the  contact  contributions  of  each  contact  node  to 
globed  tangent  stiffness  and  residual.  In  this  process,  the  change  of  the  profile  or 
bandwidth  must  be  taken  into  account. 
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The  following  vectors  are  now  introduced: 

iV^  =  {-(l-On,  -en,  n} 

T^  =  {-(l-Oei,-eci,  ei} 
t/^  =  {-n,n,0} 
du  =  {dUi,  dU2,  dUa) 
where  Tl  and  Ci  are  shown  in  Figure  B.l.  The  nodal  values  dUi,  dU2,  and  dU,  of 
the  displacement  increments  associated  with  the  three  nodes  (1,  2,  s)  are  involved 
in  the  contact  of  a  single  slave  node  with  a  master  segment.  With  the  definition  of 
these  vectors,  the  contact  residual  and  stiffness  are  given  respectively  as 

U  =  XN  (J3.25) 


(B.26) 
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Symmetry  Treatment  and  Contact  Node  Search  Strategy 


In  a  finite  element  setting,  given  the  slave  node,  the  contact  search  problem 
is  to  search  for  the  closest  projection  onto  the  master  surface.  There  are  three 
steps  to  tahe  for  the  search.  The  first  step  is  to  determine  which  master  node  is 
closest  to  the  slave  node.  In  computer  science,  this  first  step  is  often  referred  to  as 
the  “nearest  neighbor  problem”  [64].  The  second  step  is  to  determine  which  of  the 
master  segments,  having  the  master  node  as  one  of  their  vertices,  is  closest  to  the 
slave  node.  The  third  step  is  to  determine  the  closest  point  of  the  master  segment 
to  the  slave  node  in  the  biunit  domain. 

The  most  obvious  but  most  expensive  method  for  the  nearest  neighbor  problem 
is  the  global  search,  which  is  to  simply  check  the  distance  between  a  given  slave  node 
and  each  master  node,  eind  choose  the  master  node  with  the  minimum.  However, 
a  simple  local  search  method  developed  by  Benson  and  Hallquist  wzis  adopted  here 
[64].  The  idea  was  to  take  advantage  of  the  elemental  connectivity  associated  with 
the  finite  element  method.  This  is  briefly  outlined  in  Figure  C.l. 

Assume  that  at  time  the  closest  master  node  M  to  the  given  slave  node  5 
is  known.  It  is  easy  to  find  a  short  list  of  potential  new  master  nodes  for  the  next 
load  step  in  the  immediate  neighborhood  of  M  by  simply  checking  the  elemental 
connectivity.  The  list  of  candidate  nodes  is  generated  by  considering  all  of  the  nodes 
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:  Regular  master  nodes 

:  Potential  closest  points  for  next  step 

:  Current  closest  point  corresponding  to  slave  node  S 
:  The  slave  node 


Figure  C.l  Description  of  a  local  search  technique. 
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defining  the  vertices  of  all  of  the  master  segments  having  M  as  a  vertex.  The  as¬ 
sumption  behind  this  idea  is  that  during  each  time  step,  the  relative  sliding  between 
the  slave  node  and  the  neighborhood  is  not  more  than  some  accovmtable  size,  such 
as  the  element  size.  The  new  neighbor  search  is  updated  in  a  new  time  step.  The 
search  always  needs  to  find  the  closest  node  to  the  slave  node  in  its  neighborhood. 
Searching  a  small  neighborhood  is  computationally  much  cheaper  th8ui  seEU'ching 
all  the  nodes.  Generally,  at  the  initial  time  step,  a  global  search  is  performed  for 
the  initial  nearest  neighbor.  The  search  neighborhood  is  updated  appropriately  as 
the  slave  node  moves  over  the  master  surface  since  the  neighborhood  is  defined  in 
terms  of  the  last  nearest  neighbor.  The  strategy  works  well  in  most  situations. 

After  determining  the  nearest  neighbor,  the  exact  segment  onto  which  the  slave 
node  is  projected  needs  to  be  determined  .  This  step  is  called  local  search,  which 
is  very  expensive  as  it  is  generally  associated  with  nonlinear  iterations.  In  general, 
there  are  several  surface  segments  attached  to  the  nearest  neighbor  corresponding 
to  the  slave  node.  As  there  may  exist  corners,  edges,  and  complicated  curvatures 
at  the  nearest  neighbor,  the  following  strategy  is  considered. 

For  all  the  potential  segments  attached  to  the  nearest  neighbor  on  the  master 
surface,  first  calculate  the  parent  domain  parameters  (^c>  »?c)  for  all  the  segments. 
Then,  determine  one  segment  with  the  minimum  therefore  with  the  min¬ 

imum  gap  g.  However,  it  is  not  guaranteed  that  (Cc»»7c)  will  fall  into  the  biunit 
domain.  Depending  on  where  the  nearest  neighbor  is  located,  and  depending  on 
the  local  curvature,  the  potential  contact  element  types  cein  be  one  of  the  following 
three:  node-to-surface  contaw:t  element,  node-to-line  contact  element,  or  node-to- 
node  contact  element  (see  Figure  C.2).  In  the  following,  the  problem  is  discussed 
in  terms  of  the  geometric  location  of  the  nearest  neighbor  eind  the  location  of  the 
projected  node  in  terms  of  parent  domain. 


Figure  C.2  Description  of  contact  element  types. 
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In  general,  the  contact  region  for  each  delanaination  is  unknown  before  solving 
the  equilibrium  equation.  Therefore,  an  assumption  must  be  made  for  the  contact 
region  of  each  delamination.  The  best  approximation  is  to  assume  that  the  contact 
region  will  remeiin  the  same  before  and  after  deformation.  This  means  that  the 
contact  region,  which  is  known  at  the  previous  iteration,  is  going  to  be  the  same  for 
the  current  iteration.  This  assumption  allows  us  to  solve  the  equihbrium  equation 
for  the  displacement.  Following  each  solution,  the  validity  of  the  assumption  for  the 
contact  region  must  be  checked.  Some  portions  of  the  delamination  may  no  longer 
be  in  contact,  or  the  portions  that  were  not  in  contact  may  now  be  in  contact. 

The  first  step  in  the  determination  of  the  contact  region  is  to  find  the  position 
of  each  slave  node  with  respect  to  its  master  element  in  the  parent  domain,  which  is 
essential  for  both  gap  calculations  and  for  judging  what  type  of  contact  is  established 
between  the  slave  node  and  the  master  segment. 

In  order  to  calculate  the  gap,  the  coordinates  (ifc,  ^c)  of  the  point  c  on  some 
master  segment  must  be  calculated.  There  are  several  elements  attached  to  the 
nearest  neighbor.  Therefore  {^cjVc)  must  be  calculated  for  each  of  these  elements. 
The  contact  point  is  defined  on  the  meister  segment  as  the  point  closest  to  the  slave 
node  and  is  calculated  by  strictly  solving  a  minimization  problem.  The  inequality 
constraints  that  bound  the  isoparametric  coordinates  between  -1  and  +1  are  not 
explicitly  imposed.  If  the  solution  of  the  unconstrained  problem  lies  outside  the 
permissible  range,  the  node  is  usually  considered  not  in  contact  with  the  segment. 
The  function  to  be  minimized  is  denoted  J,  the  location  of  the  slave  node  is  Xc, 
and  the  vector  to  the  isoparametric  slave  segment  is  X,.  The  problem  is  stated  as 


1 

Min  :  J  =  -(*,  -  Xc){x,  -  Xc), 
J,(  =  ei(x,-Xc)  =  0, 
J,„  =  e2-ix,-Xc)  =  0 


(C.l) 
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An  explicit  solution  of  these  equations  is  very  cumbersome  emd  therefore  they  are 
solved  numerically  using  Newton’s  method.  Numerical  tests  suggest  that  in  cases 
where  elements  are  not  distorted,  two  iterations  are  sufficient  to  obtain  an  engi¬ 
neering  accuracy,  as  a  good  guess  always  exists  from  the  previous  solution  step. 
However,  the  method  diverges  with  distorted  elements  imless  the  initial  guess  is 
accurate.  A  fast,  robust  and  reasonably  accurate  approximate  contact  point  cal¬ 
culation  [64]  was  used  because  of  the  instability  of  the  Newton-Raphson  method. 
An  exact  contact  point  calculation  is  critical  in  the  delamination  contact  problem 
to  prevent  the  solution  from  oscillating  divergence  when  locating  the  contact  point. 
Two  to  three  iterations  with  a  least-squares  projection  were  used  to  generate  aa 
initial  guess  [64]. 

After  finding  the  f^he  following  tests  need  to  be  applied: 

1.  The  Contact  Gap  Test: 

The  slave  node  s  may  stay  away  from  its  master  element  (Figure  C.4.a)  if 

g  >  t,  (C.2) 

is  satisfied,  where  €g  is  the  extending  tolerance  that  is  used  to  overcome  nu¬ 
merical  errors  in  the  czdculation  of  the  gap  g  and  to  extend  the  parent  domain 
to  overcome  the  discontinuous  slope  between  elements.  In  this  case,  the  slave 
node  s  is  not  in  contact  with  its  master  element.  If  the  slave  node  s  were 
in  contact  with  its  master  element  previously,  then  the  contact  would  be  re¬ 
leased  without  applying  the  contact  force  test.  In  this  case,  the  value  of  the 
contact  force  for  the  mzister  element  of  that  slave  node  is  set  to  zero.  If  two 
delamination  surfaces  axe  not  in  contact,  there  are  no  surface  tractions  on  the 
delamination  surface. 
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% 


a.  Nearest  neighbor  is  a  comer  node. 


NO  CONTACT 


NTS  CONTACT 


NTL  (AB)  CONTACT  NTL  (CD)  CONTACT 

c.  Nearest  neighbor  is  an  interior  node. 


Figure  C.3  All  the  possible  locations  of  a  slave  node  with  respect  to  a  master 
ellement  in  the  contact  search. 
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The  slave  node  s  may  penetrate  into  its  master  element  (Figure  C.4.b)  if 

5  <  0  (C.3) 

is  satisfied.  If  the  slave  node  s  penetrates  into  its  master  element,  contact  is 
established  without  applying  the  contact  force  test. 

2.  The  Contact  Force  Test: 

The  contact  force  test  is  applied  if 

0  <  Isl  <  e,  (C.4) 

If  two  delamination  surfaces  are  in  contact,  then  the  contact  forces  must  be 
compressive.  If  the  nodal  contact  force  at  slave  node  number  S  is  foimd  to  be 
tensile,  then  the  slave  node  must  no  longer  be  in  contact  and  it  is  released.  In 
this  case,  the  value  of  the  Lagrange  multiplier  for  the  master  element  of  that 
slave  node  is  set  to  zero. 

3.  The  Geometric  Contact  Check  to  Determine  Potential  Contact  Element  Type: 
If  there  is  overlapping,  the  following  checks  must  be  made  to  evaluate  what 
type  of  contact  element  is  required: 

(a) .  Nearest  neighbor  is  a  comer  node  on  the  master  surface  (Figure  C.3. a). 
The  slave  node  s  may  project  onto  the  master  element  (Figure  C.3. a)  if 

iec|<l  +  €„  |»?c|<l  +  e,  (C.5) 

Otherwise,  no  contact  occurs  for  the  slave  node. 

(b) .  Nearest  neighbor  is  an  edge  node  on  the  master  surface  (Figure  C.3.b). 
Depending  on  the  local  coordinates  of  the  master  element,  there  are  several 
cases  to  be  considered.  For  the  coordinates  shown  in  the  figure,  if 


(C.6) 


I^cl  >  1  +  c. 
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a.  Slave  node  S  is  penetrating  into  its  master  element;  CONTACT 


b.  Slave  node  S  is  away  from  its  master  element;  NO  CONTACT 

_L 


Eg 


T 


c.  Slave  node  S  is  in  CONTACT  with  its  master  element  if  contact  force 
is  compressive 

Figure  C.4  Determination  of  contact  according  to  the  position  of  the  slave  node 
with  respect  to  its  master  element. 
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no  contact  occurs.  If 

l^cl  <  1  +  |»?cl  <  1  +  Cj  (C'.7) 

a  node-to-surface  contact  element  may  be  assumed,  depending  on  the  further 
contact  check.  If 

l^cl  <  1  +  l»7c|  >  1  +  (C.8) 

a  node-to-line  contact  element  may  be  assumed,  depending  on  the  further  con¬ 
tact  check,  with  line  AB. 

If  the  local  coordinates  reverse  in  ^  and  tj,  Eq.  (C.8)  only  needs  to  reverse 
and  Tjc- 

(c).  Nearest  neighbor  is  an  interior  node  on  the  master  surface  (Figure  C.3.c). 
There  are  four  checks  typical  for  this  case,  which  may  result  in  three  different 
contact  elements.  If 

l^cl  <1  +  65,  |t7c|  <  1  +  Cj  (C'.9) 

a  node-to-surface  contact  may  be  assumed.  It  may  result  in  three  different 
contact  elements.  If 

l^cl  >!+€,,  |t7c|  >  1  +  fg  (C.IO) 

a  node-to-node  contact  is  assumed,  with  the  closest  projection  point  being  the 
nearest  neighbor  node.  If 

|^c|>l  +  e„  |t7c|<1  +  c,  (C.ll) 

a  node-to-line  contact  may  be  assumed,  depending  on  the  further  contact  check, 
as  some  point  on  line  AB  is  the  closest  point  projection.  If 


l^cl  <  1  +  l^cl  >1  +  ^9 


(C.12) 
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a  node-to-line  contact  may  be  assiimed,  depending  on  the  further  contact  check, 
as  some  point  on  line  CD  is  the  closest  point  projection. 

After  determing  the  master  segment,  the  gap  g  for  those  node-to-line  and  node- 
to-node  contact  elements  must  be  redefined.  The  gap  emd  the  contact  force  should 
also  be  rechecked. 

The  implementation  of  the  above  treatment  is  fairly  complicated.  There  is 
another  optional  approach  to  avoid  the  complicated  treatment,  though  it  may  fail 
in  some  cases.  That  is,  by  slightly  increasing  Cg  (i.e.,  to  extend  the  parent  domain), 
node-to-surface  contact  elements  may  be  used  alone.  In  most  cases,  the  numerical 
results  differ  very  slightly  from  those  obtained  by  the  complicated  logic.  Both 
options  were  used  alternatively  in  the  simulations. 
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Contact  Residual  and  Stiffness 
of  a  3-D  Contact  Element 


In  this  section,  attention  is  restricted  to  one  contact  element,  realizing  that  the 
total  contact  contribution  is  merely  an  assembly  of  all  such  active  contact  elements. 
The  objective  is  to  derive  expressions  for  contact  residual  ^  zind  consistent  contact 
stiffness  kc-  The  contact  element  is  defined  such  that  it  contains  the  master  nodes 
of  the  segment  and  the  active  slave  node,  as  shown  in  Figure  D.l.  In  this  fig’ire, 
the  current  position  of  a  slave  node  is  given  as  and  the  master  element  surface 
is  the  element  containing  its  projection  ®c  on  the  master  surface.  As  indicated 
in  the  figure,  it  is  assumed  that  Xc  lies  in  the  continuous  interior  of  the  master 
element  siurface,  as  indicated  in  the  figure.  The  case  where  ®c  lies  on  a  discontin¬ 
uous  portion  of  the  master  surface  can  occasionally  be  troublesome.  In  practice, 
steps  may  be  tjdcen  to  avoid  this  situation  by  either  extending  element  paxameteri- 
zations  outside  the  normal  domain  (i.e.,  by  allowing  to  lie  slightly  outside  the 
bivmit  square),  or  by  altering  the  contact  element  types  near  corners,  edges,  and 
vertices  with  complicated  curvatures  such  that  the  closest  point  projection  may  be 
best  approximated.  In  the  numerical  verification  the  two  approaches  were  used 
alternatively,  with  the  former  and  less  costly  approach  being  preferred. 

First  consider  a  four  node  quadrilaterid  segment  of  arbitrary  shape.  Ftom 
this  the  three  node  triangular  version  can  be  easily  derived.  If  defining  standard 
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isop2ir2tmetric  shape  fimctions  on  the  biiinit  square  (for  3-D  problem),  the  following 
isoparametric  interpolation  of  the  geometry  of  the  master  segment  can  be  assumed 
as 

X  =  ^  =  V’o®*  (-D-l) 

a=l 

where  denotes  the  position  vector  of  a  master  node  and  (a  =  1  —  4)  are  the 
familiar  bi-linear  interpolation  functions  defined  as  follows 

=  J  (1  +  ((.)  (1  +  VI.)  (0.2) 


where  a  =  1  —  4  and  there  is  no  summation  on  a. 

Given  an  element  surface,  the  tangents  to  ^  and  ij  are  given  here  in  terms  of 
the  isoparametric  mapping: 


1 

Cl  =  x,(  =  X^  =  -ia{l  +  r]Tja)  31“ 
62  =  aJ,„  =  V’a.,  ^  »7a  (1  +  Ua)  ^ 


(D.3) 


The  normality  conditions  of  the  gap  are  as  follows: 


Ci-iXs-Xc)  =  0, 

e2-(«,-®c)  =  0 


(DA) 


Consider  now  the  total  differentials  of  the  two  equations  in  Eq.  (D.4): 


dei-(a;,-®c)  +  ei-(d«, -<iaJc)  =  0  (D.5) 

de2-{Xs-Xc)  +  e2-{dXg-dXc)  =  0  (D.6) 

Prom  Eq.  (D.3),  one  obtains  the  differentials  of  the  tangents 

de\  =  x*dr}  -}-  daJc,(  {D.7) 

de2  =  x*d(  +  dXc,,,  {D.8) 
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where  the  vector  aj*  is 

-ic^  +  a;* -aj^)  (D.9) 

which  is  interpreted  as  an  indicator  of  shape  of  the  element.  It  is  easy  to  verify 
that  for  a  plane  element  of  parallelogram  shape  aj*  =  0.  For  the  second  term  of 
Eq.  (D.5)  and  Eq.  (D.6)  one  obtains  the  following 

dx,  -  dXc  =  dU,  -  dUc  -  Cid^  -  e2drj  (D.IO) 

Finally,  introducing  the  four  equations  (Eq.  (D.7),  Eq.  (D.8),  Eq.  (D.9),  and  Eq.  (D.IO)) 
into  Eq.  (D.5)  and  Eq.  (D.6)  and  solving  for  the  increments  of  the  surface  coordi¬ 
nates  results  in 

~  ^{9{m22n  •  dUc,i  -  mijn  •  dWc.i,)  +  (m22Ci  -  >711262)  •  {du,  -  dUc)} 

(D.ll) 

dn  =  ^{^(mnn  •  dUc,n  -  •  dl*c,{)  +  (mne2  -  >711261)  •  (dit,  -  dUc)} 

(D.12) 

with  the  abbreviations 

ffii2  =  mi2  -  X*  ■  d  (D.13) 

and 

h  =  mii>7i22  -  >^12  (-0.14) 

where  the  normal  n  is  given  as  follows 


61  X  62 

||ei  X  62! 


(D.15) 


Having  defined  the  kinematical  expressions,  the  expressions  of  first  and  second  vari¬ 
ations  of  the  gap  g  must  be  derived,  which  involves  much  algebrsiic  manipulation. 
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For  simplicity,  attention  is  limited  to  the  main  equations. 

As  was  shown  before,  the  weaJc  form  of  the  equilibrium  equation  requires  the 
first  variation  of  the  gap  g  (Eq.  (3.9))  with  respect  to  the  node  variables.  Taking 
the  first  variation  of  g  results  in 


6g  =  n-(6u,-6Uc)  +  Sn-{x,-Xc)  (D.16) 

As  6ti  is  perpendicular  to  n,  the  second  term  of  6g  vanishes.  The  second  variation 
of  the  gap  is  then 


d{6g)  =  dn  •  -  8Uc)  +  n  •  d{6u»  -  6Uc)  (-0.17) 

Some  algebraic  manipulations  show  that  the  variation  of  71  results  in 

dn  =  ^(/- n(8)n)d(ei  X  62)  (JD.18) 

I  is  the  identity  matrix  eind  0  denotes  tensor  product.  Q  is  the  surface  element 
area  and  is  defined  as 


ft  =  ||eixe2||  (£>.19) 

It  is  obvious  that  we  need  to  derive  d(ei  x  62).  By  using  Eq.  (D.7)  and  Eq.  (D.8), 
one  obtains 

d{ei  X  62)  =  {X*  X  e2)dT]  -  («*  X  ei)d^  +  {du^,^  x  62)  -  {dUa,,,  x  eO  (P.20) 
where  dr}  and  d^  can  be  eliminated  by  applying  equations  Eq.  (D.ll)  and  Eq.  (D.12). 
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After  some  algebraic  manipulations  the  variation  of  n  becomes 


dn  =  — (/  -  n  0  n){(dUc,{  x  C2  +  ei  x  dUc,^) 

+  i[(ei  X  X*)  0  (1712261  -  mi2e2)(<i«,  -  dUc) 
+  (x*  X  C2)  0  (7711162  -  mi2ei){du,  -  dUc)] 

+  ^[(61  X  ®*)  ®  n(m22dUc,^  -  mi2dUc,r,) 


(JD.21) 


+  (x*  X  62)  0  n(miidUc,,  -  mndUc,^)]) 

It  remains  to  evaluate  the  second  term  of  Eq.  (D.16),  which  is  simply  given  by 


d(6u^  -  SUc)  =  -{6Uc,idi  +  SUc,j,dr})  (D.22) 

By  applying  Eq.  (D.ll)  and  Eq.  (D.12)  to  eliminate  d^  and  drj,  one  obtains 

1 

n  •  d(6Ua  -  SUc)  =  -■^[<5Uc,{(7n22(”'  ®  61)  -  mi2{n  0  e2))(d'U,  -  dUc) 

+  SUc,r,{mu{n  0  62)  -  7ni2(n  0  ei))(du,  -  dUc)] 

-  |[^Uc.«(n  ®  n){miidUc,^  -  mi2dMc,,) 

+  SUc,r,{n  0  n)(7nii<iMc,,  -  fhi2dUc,()] 

{D.2Z) 

The  matrix  for  contact  contribution  is  expected  to  be  symmetric.  However,  only 
the  last  part  of  Eq.  (D.23)  shows  symmetry  a  priori.  To  obtain  the  symmetric 
expressions  requires  additional  algebraic  operations,  summarized  in  the  following. 
First  the  contribution  due  to  warping  must  be  separated  from  the  other  peirts  by 
considering  X*  with  respect  to  the  local  tangential  system  61,  62  ,  and  n  according 
to 

a;*  =x*ei +x^e2  +  x’n  (D.24) 

Then  it  is  necessary  to  evaluate  the  cross  products  which  give  the  corresponding 
contravariant  base  vectors.  In  each  case,  a  set  of  dual  basis  vectors  may  be  computed 
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such  that  e'  •  Cj  =  ^ j  (i,  j=l-2).  A  dual  basis  at  the  contact  point  Xc  may  be 
similarly  defined  as 

Ca .  =  si  (D.25) 

A  merit  tensor  may  also  be  expressed  as 

=  Ca  •  (£>.26) 

The  contravariant  tangents  may  be  defined  as 

e“  =  (ma^)“'e^  (£>.27) 


After  rearranging  the  terms  in  a  convenient  manner,  the  final  equation  is 


x^Q 

d(Sg)  =  (-^mi2  -  l)[dite,{(n  0  6^)  +  dUc,„{n  0  e^)](6u,  -  SUc) 
+  -  l)[^Uc,{(n  0  e^)  +  6Uc,^(n  0  e^)](du,  -  dUc) 

-  0  n)SUc,(  +  miidUc,n{n  <S>  n)SUc,,,)] 

-  7fii2(SUc,((n  0  n)dUc,,,  +  SUc,r,{n  0  n)dUc,^] 

3.3 

-  -  dUc)[{miim22  -  mj2mi2)((e^  0  e^)  +  (e^  0  6^)) 

+  i^y(m22(e^  0  e^)  +  mii(e^  0  e^))]{6u,  -  SUc) 

Q 

— ^[m22dUc,{(n  0  e^)  +  miidUc,,,{n  0  e^)](6u,  -  SUc) 
n 

X^Q 

-  -^[m22SUc,^{n  0  e^)  +  mn6Uc,,(n  0  e^)](du,  -  du^) 


(£>.28) 


which  proves  the  symmetry.  It  is  desired  to  replace  the  displacement  vectors  by 
introducing  the  displacement  interpolation  according  to  Eq.  (D.l).  The  second 
order  part  of  the  contact  stiffness  then  follows  from  the  identity 


d{6g)  =  du^ 


d'^g 

du^du^ 


6u^ 


(£>.29) 
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From  the  virtual  work  linearization,  the  contact  stif&iess  can  be  written  as 


1  dg  dg  d  ,  dg  . 


(D.30) 


And  the  contact  residual  force  of  the  element  is 


f;  = 


(11.31) 


For  computation2d  purposes  it  is  advantageous  to  rewrite  the  result  in  matrix  no¬ 
tation.  To  achieve  this,  it  is  first  necessary  to  mrdce  the  following  definitions  for 
vectors  such  that  the  slave  node  is  represented  as  the  last  node: 


=  {-«  >1,  -nV’a,  •  • ' ,  n} 

V/  =  {-eVi,  -6^2,  •  •  • ,  e^} 
=  {nV'i.f,  n^2,e,  •  •  • ,  0, 0, 0} 
UJ  =  {nV'i.v,  nxi?2,r,,  •  0, 0, 0} 

du  =  {dui ,  dti2,  •  ” ,  du,} 


(£>.32) 


It  is  noted  that  all  the  vectors  are  bjised  on  the  interpolation  functions  evaluated  for 
the  surface  coordinates  of  point  c.  Using  these  vectors,  the  contact  reaction  forces 
(contact  residual  forces)  on  the  contact  element  with  an  active  slave  node  are  then 
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The  corresponding  consistent  tangent  stiffness  matrix  takes  the  form 

fee  =  €NN^  +  tK^rhn  -  l)[t;.  v,’-  +  V,V/  +  V. [//■  +  V,un 

-  -  m,2(U,U^  +  t/jt/f)) 

-  V/  +  Vj  V^)  +  V,  +  m^VjV/)) 

-  V/  +  Vat/f)  +  V,^  +  V.I/f )]} 

ft 

(D.34) 

where  the  t  =  eg  for  the  penalty  method,  and  <  =  A  +  egf  for  augmented  Lagrangian 
treatment.  It  must  be  emphasized  that,  except  for  the  linear  term  NN^,  all 
terms  are  due  to  the  second  variation  of  the  gap  and  therefore  accoimt  for  its 
nonlinear  kinematics.  It  is  seen  that  the  contribution  due  to  warping  represented 
by  those  parts  containing  is  well  sepzirated  from  the  others.  Obviously,  for  linear 
kinematics,  warping  does  not  enter  into  consideration. 

As  a  special  case,  the  residual  and  stiffness  matrix  for  a  three  node  triangular 
contact  segment  can  be  extracted  by  simply  letting  x®  =  0  in  Eq.  (D.34).  Another 
simplification  is  due  to  the  fact  that  the  physical  point  of  contact  can  be  explicitly 
calculated  since  the  normal  n  is  constant.  Expressions  for  node-to-line  element  and 
node-to-node  element  in  3-D  space  can  be  similEurly  derived.  Expressions  for  planar 
and  axisymmetric  problems  ran  be  readily  derived  using  an  approach  similar  to  the 
one  used  above. 

For  the  case  where  one  body  is  deformable  and  the  other  is  rigid,  and  the  motion 
of  the  rigid  body  (master  surface)  is  prescribed,  it  is  assumed  that  the  surface  can 
be  represented  either  by  a  set  of  element  edges  defining  the  surface  or  in  a  globally 
continuous  fashion.  The  effect  of  this  assumption  is  simply  to  eliminate  the  master 
degrees  of  freedom  from  the  problem,  leaving  only  the  slave  degrees  of  freedom  as 
unknowns. 
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Input  Data  Files  for  Computer  Codes 


§E.l  Interface  for  2CRACK 

The  model  described  in  Chapter  3,  was  implemented  into  a  nonlineeu:  finite 
element  program  designated  as  “2CRACK”.  An  user-fiiendly  interface  to  2CRACK 
was  written  to  enable  users  unfamiliar  with  finite  element  procedures  to  analyze 
laminated  composite  panels  either  with  or  without  pre-existing  matrix  cracks  and 
delaminations. 

A  typical  terminal  session  using  2CRACK  is  listed  below.  The  problem  solved 
is  a  flat  panel  subjected  to  a  line  load  as  described  in  Chapter  3. 
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g  _ »g‘6o  g66f-onv-og _ T-"j.Nrj/)ciNilaNvana|  TvnuSMvxs  I  i  ^g  f»0  g(.Gi  uiv  og  ! '.wi  xn.wi  lUNVMfuii  ivj>i$hvxs 


U>AD1N6  ••••••«••**••*•••••*»•••  1  Mould  you  like  to  read  sone  information  about  boundary  conditions?  ’y/n) 


STAR$rYIAl  :  inURANDlIM'UT  im'.  l  20- AUG -1992  09:24  P^ge  S  LSTAR$DUA1 :  (DURANDj  INPUT.  IlfT,- 1  20-AUG-1992  09:24  Page  6 
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10  fX)NTINUE  PRESS  •‘KB’niKN" 


STAM$1)IIA1 :  |WIRANl)|  lNI>lfr- ir/r.  1  20  AUU  1902  09:24  Pdqe  7  I  I  STAR$DOAl :  (DURANDI  INPUT.  IW;  1  20-AUG-1992  09:24  Page  8 
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,STAR$DUA1: (DURANDIINPUT.DAT;!  20-AUG-1992  09:24  Page 


fl  mesh:  flat:  “fl",  cylindrical  :  "cy" 

I  0  ncrack: 1-originally  cracked;  0-originally  perfect 

[-0. 40000E-02  du:  displacement  increment 

j  0.15000E+04  xload:  maximum  load 

|150  itter:  load  step  number  limit 

0.15200E+01  arc:  arch-arc  angle .plate-model  length 

O.OOOOOE+00  rin:  inner  radius:  plate:  0.0 

1  0.89600E-01  hi;  plate  thickness 

!  O.lOOOOE+00  hO:  plate  thickness  with  strain  gage  attached 
I  3  ibound:  CC.R:3;  MSS.R;9;  CC.RN:10,  SS.RN;11.  etc 

'  1  nip:  SS  point  location  for  IB0UNI>9 

i  0.17600E+08  E_xx 

i  Q.14100E+07  E_yy 

;  0. 14100E+07  E_Z2 

'  0.81000E+06  G_xy  for  linear  case,-  G*o_xy  fomonlinear  case 

■'1  50360E+06  G_y2 

O.SlOOOE+06  G_x2 

0.29000E+00  Nu_xy 

0.40000E+00  Nu_yz 

!  0.29000E+00  Nu_xz 

0  In-situ  strength?  l:yes;  O.-no 

2  icrit : 1-Tsai;2-Hashin, 3-Gosse;4-Chang-Choi 

I  0.22000E+06  X_t 

I  0.23100E+06  X_c 

'  0.64500E+04  Y_t  for  constant;  Y*o_t  for  in-situ 

C.13000E+01  A  an  empirical  constant  for  in-situ  Y_t 

0.70000E*00  8  an  empirical  constant  for  in-situ  Y  t 

0.36700E+05  Y_o 

C,155C0E+O5  Sliear  strength 

0 . 2r'000E*01  C  an  empirical  con.stant  for  in-situ  S_c 

0 , 1000CE*01  D  an  empirical  constant  for  in-situ  S  c 

0.0C00OE*0O  Alpha 

0.50PO0E+0O  Strain  energy  release  rate  for  mode  I,  G_Ic 

,  0.i8000E+01  Strain  energy  release  rate  for  mode  11,  G_IIc 

3  Total  number  of  stacks 

3  Number  of  plies  in  stack  no.  1 

O.OCOOCE+GO  Orientation  of  stack  no.  1 

3  Number  of  plies  in  stack  no.  2 

O.90G00E+O2  Orientation  of  stack  no.  2 

j  3  Number  of  plies  in  stack  no.  3 

;  O.OOOOOE+00  Orientation  of  stack  no.  3 


18 


i  4 

6 

' loaddis . dat 
I'  drivel  .dat 
i'drlve2  .dat 
i'  stnl .  dat 
' stn2.dat 
''mesh,  dat 
'  fail . dat 


naply:  Total  number  of  actual  plies 

nelarc:  number  of  elements  on  the  length  (or  circumference) 
nmdat:  number  of  increments  interval  you  want  your  output 
igrow:  1: initial  damage  growth;  0:no  growth 
nlkrall:  location  of  matrix  crack  (lower) 
n3kral2:  location  of  matrix  crack  (upper) 
output  LOAD-DISPLACE.MENT  FILE  e. g. loaddis. dat 
plotting  data  file  name  FOR  #1  DELAMINATION, e.g.  drivel.dat 

plotting  data  file  name  FOR  #2  DELAMINATION.e.g.  drive2.dat 

plotting  data  file  name  FOR  #1  DELAMINATION.e.g.  stnl.dat 

plotting  data  file  name  FOR  #2  DELAMINATION.e.g.  stn2.dat 

deformed  shape  data  file  name,  e.g.  mesh.dat 
FAILURE  data  file  name,  e.g.  fail.dat 


1 

I 
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The  model  described  in  Chapter  4,  was  implemented  into  a  nonlinear  finite 
element  program  designated  as  “3CRACK”.  An  input  data  and  the  Uesr  Manual 
to  3CRACK  were  written  to  enable  users  unfamiliar  with  finite  element  procedures 
to  analyze  laminated  composite  panels  either  with  or  without  pre-existing  matrix 
cracks  and  delaminations. 

A  typical  input  data  file  using  3CRACK  is  listed  below.  The  problem  solved 
is  a  fiat  panel  containing  an  existing  surface  crack,  Ein  internal  shear  crack  and  an 
elliptical  delamination  subjected  to  a  sphericed  indenter  as  described  in  Chapter  4. 
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SAMPLE  INPUT.DAT  DELAMINATION 
WITH  BOTH  TYPES  OF  CRACKS 


_STAR$DUA1: [DURAND] INPUT. DAT; 2  20-AUG-1992  09:35  Pagel 

0  igrad;  consider  updating  stiffness  1;  yes  0:  no 

0  liucd:  consider  material  degradation  1:  yes  0:  no 

-2  load:  load  type  0:  POINT;  1:I<INE;  2:CONTACT  3:Beam  Comp;  -2:  Hertz  cont 

-0.005  du:  DISPLACEKENT  INCREMENT 

0.25  rimpact:  indenter  radius 

250.  xload;  MAXIMUM  lOAD  desired 

3  IGEO:  Geometry  type  1: CYLINDER;  2. ARC;  3: BEAM 

9  IGEOS:plate  type  l:solid;  2:DCB;  3:  plate  with  ec  del;4:  del+shear  cr 
1  MOPTN:l:  1/4  model  {1/2  for  DCB) ;  2:1/2  model  (full  for  DCB) 

2.0,1.5,0.104  xleng,yleng,zleng:  sizes  of  the  plate 

0.0  Rin:  inner  radius  for  the  curved  plate,  zero  for  flat  panels 

+18  IBOUND:l:CC.FF  2:CC.SYM  3 : CORNER . CC . SYM  7.BE.CC.SYM  18:CC+bC+SC  19:  20: 

9  nxsegt:  total  H  of  secjnents  in  x  direction 

1.1. 1.1. 1.1. 2. 1.1  nxseg(l-nxsegt) :  I  of  elements  in  each  segment. 

!o. 03, 0.06, 0.02, 0.01, 0.01, 0.01, 1.19, 0.2, 0.47  xseg(l-nxsegt) : segments  sizes 
14  nssegt;  total  »  of  segments  in  theta (s)  direction 

ll,  1.1,1  nsseg(l-nxsegt) :  »  of  elements  in  each  segnent. 

22.5.22.5.22.5.22.5  sseg(l-nssegt) : segments  in  degrees 

4  nzsegt:  total  »  of  segments  in  thickness  direction 

a,  1,1,1  nzseg(l-nxsegt) :  »  of  elements  in  each  se^nent. 

0.013,0.039,0.039,0.013  zseg(l-n2segt) :  segments  sizes 

I  NDEL;  number  of  delaminations,  limited  to  1,  or  0  for  undelaminated  panels 

5.1  nxseoc.nzsegc;  delamination  location  *  of  segments  in  x  and  z 

C. 10,0.1510.0899  ra,rb,rc:si2e  of  del . long/short  ax,rc<min(ra,rb)-nxseg(nxsegc) 
0.013  dzlc:  delamination  location  in  thickness  direction 
,2,0.375,0.013  icshape(l)  ,rxc(l) ,r2c(l)  :bend  crack,  type,  x  and  z  sizes 
0.151  ryc(l):  shear  crack  length  (in  y) 

1.0.03  nbshcr.shcr:  location  in  segment  »  and  size  of  shear  crack  in  x 
0.0131,0.0905  shzcl,shzc2;  location  of  shear  crack  in  z 

i23.2e+06  ,  1.3284e+06  ,  1.3204e+O6  ex,ey,ez:  young’s  modului  in  x,y,and  z 

0.9096e+06  ,  0.9086e+06  ,  0.90a6e+06  gxy , gy z , gxz :  shear  mogului 

!.29  ,  .28  ,  .28  vxy,vyz,vxz:  poisson’s  ratios 

>412930.  ,190840.  XX,  XXC:  tensile  and  compressive  strength  in  x 

6410., 24346.  YY,  YYC:  tensile  and  compressive  strengtJi  in  y 

,5877.  SSO:  in-plane  shear  for 

;0.0  ALPHA:  SHEAR  NONLINEARITY  FACTOR 

jo.,0.  AXX,  AYY:thermal  expansion  coefficients  in  x  and  y 

;75.  TR;  Stress  free  temperature 

|75.  TO:  Room  temperature 

>2  icrit:  type  of  criterion  for  initial  damage:  l:Tsai  2;Hashin  3:Gosse 

1.5  ,18.0,18.0  Gmaxn , gmaxs , gmaxt :  fracture  toughnss  in  Mode  I, II, and  III 

k  NPLYG:#  OF  PLY  GROUPS,  refer  to  nzsegt,  counting  from  the  bottom 

|2  «  OF  LAYERS  IN  »1  GROUP 

;0.  PLY  ORIENTATION  OF  «1  GROUP 

II  »  of  layers  in  group  #2 

90.0  ply  orientation  of  ply  group  *2 

1  #  of  layers  in  group  13 

90.  ply  orientation  of  ply  group  #3 

1  »  of  layers  in  group  14 

0.  ply  orientation  of  ply  group  #4 

1.1.2.1.1.1.1.10000.  l,icontac,islt,npas,isech,ictype,n4ug,penp 

2.1.2.2.1.1.1.10000.  2,icontac,islt,npas,isech,ictype,naug,penp 

io. 
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SUMMARY 

This  manual  describes  the  input  data  used  by  3CRACK.  For  convenience,  the 
input  data  is  read  from  one  hie:  input.dat.  The  input  data  is  arranged  in  the 
following  sequence: 

1)  Global  parameters  and  control  data. 

2)  Loading. 

3)  Geometry  of  the  panel. 

4)  Boundary  conditions. 

5)  Finite  element  mesh. 

6)  For  cracked  specimens,  the  location  of  the  cracks. 

7)  Material  properties. 

8)  Lay-up. 

9)  For  contact,  specify  contact  stiffness  and  solution  algorithm. 

The  data  is  read  in  format  free  input;  hence,  comments  can  be  inserted  in  the 
input  hies  between  the  data  lines  and  following  the  last  item  in  a  data  line.  Here,  a 
data  line  is  dehned  as  a  sequence  of  input  values  that  follow  continuously,  usually  on 
one  line. 
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n^IPUT.DAT-a  manual 

i.  Control  Parameters 

igrad:  updated  stiffness  index 

1:  the  material  stiffness  coefficients  Cij  are  updated.  In  general,  the  convergence 
can  be  improved; 

0:  Cij  =  Cij  is  assumed.  It  can  save  computing  time  in  calculting  the  Cij  due  to  the 
deformation. 

lined:  material  degradation  index 

1:  Material  degradation  model  is  activated; 

0:  No  material  degradation  is  considered.  It  is  suggested  not  to  consider  material 
degradation  for  the  time  being. 

£.  Loading  Parameters 
load:  load  type  index. 

0:  a  point  load  applied  at  the  center  of  the  plate; 

1 :  a  line  load  uniformly  applied  in  the  center  line  of  width  direction; 

2:  contact  load  applied  by  in-situ  cont^u:t; 

3:  axial  compression  applied; 

-2:  Hertzian  contact  applied. 

dn:  displacement  load  increment  for  each  step. 

ndeb:  line  load  type  index.  Not  used  for  other  types  of  load. 

2:  a  line  load  applied  on  upper  and  lower  center  line  of  the  specimen  such  as  in 
DCB  I  loading  ; 

1:  a  line  load  applied  on  upper  center  line  of  the  specimen  such  as  in  DCB  II,  I+II, 
or  cylindrical  line  loading  ; 

-1:  a  line  load  applied  on  the  lower  center  line  of  the  specimen. 
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rimpact:  indenter /impact  radius.  Enter  only  for  contact  loads  (Load=2,  or  -2). 
xload:  the  maximum  load  desired. 

S-Geometry  Parametera 
igto:  geometry  type. 

1 :  for  a  full  cylinder; 

2:  for  a  arch  ; 

3:  for  a  beam  or  plate. 

igeos:  plate  type. 

1 :  for  a  solid  beam  or  plate; 

2:  for  a  DCB  beam  (with  one  through-the-width  crack); 

3:  for  a  plate  with  an  elliptical  delamination  (with  or  without  a  bendingcrack); 

9:  for  a  plate  with  an  elliptical  delamination  offset  by  an  internal  shear  crack  (with 
or  without  a  bending  crack). 

moptn:  symmetry  index  for  analysis. 

1:  1/4  model  for  regular  plate,  1/2  for  a  DCB  or  a  cantilever  plate; 

2:  1/2  model  for  a  reqular  plate,  full  model  for  a  DCB  or  a  cantilever  plate.  In 
most  cases,  half  or  1/4  model  is  used. 

xleng,yleng,zlcng:  beam  or  plate  dimensions  in  the  length,  width,  and  thickness 
directions. 

arc:  arch  angle  in  terms  of  degrees  (enter  only  for  arc  and  cylinders.) 
rtn:  for  arc,  enter  the  inner  radius;  for  plate  enter  0. 

4- Boundary  Conditiona  (aupport  conditiona) 

ihound:  boundary  conditions  type  (including  matrix  cracks  conditions) 

1:  left  side-full  clamped,  other  three  sides-free  (typical  for  cantilever  beam  of  DCB); 
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2:  back  and  front  sides  -  free,  left  side  -  full  clamp,  right  side-symmetry  roller 
support; 

3:  back  and  front  sides  -  free,  left  side  -  upper  and  lower  edges  compressive  clamp 
(allowing  the  axial  movement),  right  side  •  symmetry  roller  support; 

7:  back  and  front  sides  -  free,  left  side  -  full  compressive  clamp  (allowing  axial 
movement),  right  side  -  symmetry  roller  support; 

16;  back  and  front  sides  -  free,  left  side  -  full  clamp,  right  side  -  symmetry  roller 
support,  bending  crack  in  the  x-z  plane  (simulating  splitting  the  outermost  0 
degree  layers); 

18;  back  and  front  sides  -  free,  left  side  -  full  clamp,  right  side  -  symmetry  roller 
support,  bending  crack  in  the  x-z  plane  (simulating  splitting  the  outermost  0 
degree  layers).  In  addition,  a  shear  crack  at  a  distance  from  the  indenter  is 
simulated; 

19;  1/4  model  for  a  four  sided  simply  supported  plate.  Back  -SS,  front  -  symmetry 
roller,  left  side  -  SS,  right  side  -  symmetry  roller  support; 

20;  1/4  model  for  a  four  sided  simply  supported  plate  but  constrained  to  have  u)  =  0 
on  the  bottom  of  the  plate. 

5’Meah  GencTation 

First  enter  the  information  needed  for  x  direction; 
nxatgt  number  of  the  segments  in  x  direction  for  mesh  generation. 
(nzstg(i),i=l,nxsegi):  number  of  the  elements  in  each  segment. 
(x3eg(i),i=l,nx3egt):  size  of  each  segment. 

(It  is  noted  that  a  non-uniform  mesh  can  be  formed,  which  is  helpful  for  areas 
such  as  concentrated  loading  zone  and  the  crack  front,  where  a  very  fine  mesh  is 
desired.  There  are  several  checks  to  be  done  to  make  sure  the  mesh  is  right.  The 
key  checks  are  the  length  of  the  beam  (xleng),  the  critical  shear  crack  location,  the 
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delamination  size.  For  example,  the  sum  of  xseg(l  to  nxsegt)  has  to  be  the  length  of 
the  plate.  More  will  be  discussed  in  the  following.) 

Then  enter  the  information  needed  for  y  or  theta  direction: 

For  igeos  >  3  and  igeos  <  9,  ent«  information  for  panels  contaiining  elliptical 
delaminations.  In  this  case,  the  sseg(i)  has  to  input  the  angles  in  terms  of  degrees. 

nssegt  number  of  the  segments  in  angle  direction  for  mesh  generation. 

(nastg(i),i=l,nsstgt)-.  number  of  the  elements  in  each  segment. 

(sseg(i),i=l,nsaegi):  size  of  each  segment,  the  sum  of  them  must  be  90  degrees. 

Otherwise,  enter  information  for  panels  without  delamination. 

nysegt  number  of  the  segments  in  y  direction  for  mesh  generation. 

(nystg(i),i-l,nyatgi):  number  of  the  elements  in  each  segment. 

(yaeg(i),i-l, nysegt):  size  of  each  segment. 

6- Delamination  and  Matrix  Cracka 

ndel:  input  delamination  index. 

0:  no  delamintion; 

1:  one  delamination; 

If  there  is  a  delamintion  {ndel  ^  0),  determine  the  location  of  the  delamination; 
nxaegc,nzaegc:  delamination  location  index. 

nxaegc:  Niunber  of  segments  in  the  span  length  before  the  delamintion  front  in  terms  of 
nzaegt  and  nxaeg(i)\ 

nzaegc:  Number  of  segments  in  the  thickness  direction  before  the  delamintion  front  in 
terms  of  nzaegt  and  nz3eg(i). 

ra,rb,rc:  delamination  size  index. 

ra:  elliptical  size  in  x  20cis; 

rb:  elliptical  size  in  y  axis; 
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rc:  a  number  used  for  modified  elliptical  transformation  in  generating  the  mesh. 
And  it  should  be  a  number  smaUer  than  either  ra  or  rb  such  that  the  one  or 
two  row  elements  in  the  'delamination  front  is  aJlowed.  The  estimated  formula  is 
given  as:  rc  <  min(ra,rb)  -  xseg(nxsegc). 

dzlc:  delamination  location  in  thickness  direction 

If  there  is  any  bending  crack  {ibound  >  14  and  ihound  <  18),  the  location,  shape, 
and  size  of  the  bending  crack  need  to  be  determined. 

icshape(l ),rxc(l ),Tzc(  1 ):  bending  crack  location,  shape,  and  size. 
ic3kape(  1 ):  bending  crack  shap:  2-rectangular;  1-elliptical.  2  is  recommened. 
rzc(  1 ):  bending  crack  length  in  x  direction. 
rzc(  1 ):  bending  crack  depth  in  thickness  direction. 

If  there  is  any  shear  crack  {ibound  =  9),  the  location,  shape,  and  size  of  the 
shear  crack  needs  to  be  determined. 

ryc(  1 ):  shear  crack  length  in  y  direction 

nb3hcr,3hcr.  shear  crack  location  in  x  direction 

nbahcr.  number  of  segment  in  terms  of  nx9tg(i). 

shcr:  location  in  terms  of  size  in  x,  which  is  the  distance  from  the  indenter  or  the 
center  line  of  the  plate. 

shzclfShzci:  shear  crack  location  in  z  direction. 
ahzcl:  the  location  of  the  lower  tip  of  the  shear  crack. 
shzcg:  the  location  of  the  upper  tip  of  the  shear  crack. 

The  key  check  here  is  that  the  mentioned  quantities  are  generedly  associated  with 
the  0/90  or  90/0  interface  heights. 

7- Material  Properties 

EX,EY,EZ:  enter  Young’s  moduli  EX,EY,EZ. 

GXY,GYZ,GXZ:  enter  the  shear  moduli  GXY,GYZ,GXZ. 

VXY,VYZ,VXZ:  enter  the  Poission’s  ratios:  VXY,VYZ,VXZ. 

XX,XXC:  enter  the  strength  in  longitudianal  direction  XX,XXC. 
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YYfYYC'.  enter  the  strength  In  transverse  direction  YY,YYC. 

550;  enter  the  shear  Strength  SS. 

ALPHA:  enter  the  nonlinear  shear  parameter  :  ALPHA. 

AXX,AYY:  enter  the  thermal  coefficients  (AXX,AYY) 

TR:  enter  the  stress  free  temperature. 

TO:  enter  the  room  temperature. 

ICRIT:entet  failure  criteria  options.l:TSAI  2;HASHIN  3:GOSSE. 

GMAXN,  GMAXS,  GMAXT:  enter  the  fracture  toughness  in  Mode  I,  II,  and  III. 

8  -  LAY-UPS 

NPLYG:  enter  the  number  of  ply  group. 

The  following  is  for  each  layer  counting  from  the  bottom  ply. 

NLYG(I):  enter  the  number  of  plies  for  ply  group  I 
ALYG(I):  enter  the  ply  orientation  for  ply  group  I 

9  •  contact  algorithm  input 

i,icontac(i),islt(i),npas(i),iaeck(i),ictype(i),  nang(i),penp(i):  contact  interface 

data 

i:  the  i-th  contact  surface.  In  general,  two  contact  surfaces  are  considered.  One  is 
the  contact  between  the  indenter  and  the  plate.  The  other  is  the  delamination 
contact. 

icontac(i):  1:  activating  the  contact  algorithm;  0:  no  contact  considered. 
islt(i):  1;  penalty  method;  2:  augmented  Lagrangian  method. 
npas(i):  1:  one-pass  algorithm;  2:  two- pass  algorithm. 

isech(i):  1:  the  type-1  simple  search  algorithm;  2:  more  rigious  search  algorithm. 
ictype(i):  contact  element  type.  0:  node-to-node;  1:  node-surface;  2:  adaptive. 
naug(i):  enter  1  please. 
penp(i):  penalty  parameter. 
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